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Introduction 


A  major  issue  facing  the  Army  is  the  restoration  of  explosive  and  energetics 
contaminated  soil  and  groundwater.  Explosives  contamination  may  result  from  a 
variety  of  activities  including:  explosives  and  energetics  manufacturing; 
munitions  maintenance,  load  and  pack  facilities;  and  live  fire  training  activities. 
Of  particular  concern  to  the  Army  is  the  wide  distribution  of  relatively  low 
concentrations  of  explosives  and  energetics  on  many  live  fire  training  and  testing 
ranges.  The  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC)  is 
actively  pursuing  development  of  low-cost  remediation  technologies  to  restore 
these  lands  and  groundwaters  to  environmentally  acceptable  conditions. 

Base-induced  transformation  (BIT)  of  explosives  has  shown  promise  as  a 
rapid,  low-cost  technology  for  remediating  explosives  in  soil  and  water.  When 
used  as  a  restoration  technology,  the  BIT  process  incorporates  the  addition  of 
highly  basic  material  to  the  soil  or  groundwater,  resulting  in  the  transformation  of 
the  parent  compound.  Ideally,  the  process  would  drive  the  degradation  of  the 
parent  to  environmentally  benign  compounds. 

Although  the  BIT  process  has  been  known  for  decades,  the  majority  of  work 
reported  on  the  phenomena  has  been  related  to  treating  much  higher 
concentrations  than  generally  found  on  sites  requiring  remediation.  The 
mechanism  and  products  of  the  Hydroxide-  (OH-)  2,4,6-Trinitrotoluene  (TNT) 
reaction  at  concentrations  of  environmental  concern  are  poorly  understood. 
Previous  kinetics  studies  of  the  OH-TNT  reaction  have  only  considered  the 
reduction  of  the  TNT  analyte  concentration,  not  the  reaction  of  intermediates  or 
final  products.  During  these  previous  studies,  TNT  quickly  degraded  in  less  than 
40  min  at  room  temperature  in  basic  solution,  but  the  overall  OH-TNT  reaction 
may  not  have  been  completed.  Early  toxicological  studies  have  shown  that  the 
final  products  of  the  OH-TNT  reaction  may  be  benign  (Hansen  et  al.  2001),  but 
intermediate  products  may  be  toxic  to  inherent  bacteria  or  bind  with 
environmental  matrices  (Michels  and  Gottschalk  1994;  Pennington  et  al.  1995). 

Identification  and  quantification  of  the  reaction  mechanism  and  all  reaction 
components  (reactants,  intermediates,  and  final  products)  of  the  OH-TNT 
reaction  are  necessary  to  determine  the  technical  feasibility  of  BIT  as  a 
remediation  technology.  The  rate  constants  for  each  step  in  the  reaction  derived 
from  the  reaction  mechanism  could  be  used  to  guide  engineers  and  technicians 
responsible  for  designing  and  operating  the  technology. 
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Detailed  study  of  the  mechanism  of  degradation  and  the  interactions  of  the 
individual  degradation  products  requires  an  analytical  technique  that  enables 
separation  of  the  individual  compounds.  Typically,  the  individual  reaction 
components  are  separated  using  techniques  such  as  liquid  chromatography  and 
identified  using  mass  spectroscopy.  Unfortunately,  the  final  products  of  the  OH- 
TNT  reaction  make  this  process  problematic  because  of  limitations  on  standard 
separation  techniques.  The  polymers  (Felt,  Larson,  and  Hansen  2001b)  that  form 
as  a  result  of  the  OH-TNT  reaction  often  foul  the  separation  column  because  of 
their  large  molecular  size.  To  overcome  this  difficulty,  the  study  reported  here 
evaluates  an  alternate  technique,  ultraviolet-visible  (UV/VIS)  spectral  analysis,  to 
obtain  similar  information  using  factor  analysis  (FA)  of  kinetic  spectral  data.  The 
knowledge  of  basic  properties  of  the  OH-TNT  reaction  presented  in  this  report 
adds  to  our  understanding  of  the  BIT  process. 
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2  Literature  Review 


Contamination  of  groundwater,  surface  water,  and  soil  by  explosives  has 
occurred  at  military  sites  throughout  the  world  as  a  result  of  manufacture  of 
explosive  compounds,  assembly  of  munitions,  and  deployment  of  explosives 
containing  devices  (Roberts  and  Hartley  1992;  Jenkins,  Thome,  and  Walsh  1994; 
Pennington  et  al.  1995).  TNT  has  been  listed  as  a  priority  pollutant,  a  Class  C 
carcinogen,  and  has  adverse  effects  on  humans,  plants,  and  animals  (Rosenblatt 
et  al.  1991;  Smith  1991;  Won,  DiSalvo,  andNg  1976;  Palazzo  and  Leggett  1986; 
Simini  et  al.  1995).  TNT  remediation  methods,  such  as  incineration,  composting, 
bioremediation,  and  photolysis,  have  been  used  with  mixed  success  (Bruns-Nagel 
et  al.  1998;  Dillert  et  al.  1995;  Lang  et  al.  1998;  Haselhorst  1999;  Hawari  et  al. 
2000).  Existing  technologies  are  hindered  by  problems  such  as  high  energy  costs, 
costly  equipment,  extensive  soil  excavation,  very  slow  degradation  rates,  or 
incomplete  degradation  of  the  explosive  contaminants. 

The  transformation  of  TNT  in  basic  solutions  has  long  been  established 
(Janowsky  1891)  and  could  potentially  be  part  of  a  rapid  and  low-cost  tech¬ 
nology  to  remediate  TNT  contamination.  The  rate  of  TNT  transformation  when 
exposed  to  ultraviolet  light  (Dillert  et  al.  1995),  to  iron  (II)  (Brannon,  Prince,  and 
Hayes  1 998),  or  to  a  combination  of  UV  light,  ozone,  and  electrohydraulic 
discharge  (Lang  et  al.  1998)  is  enhanced  at  alkaline  pH.  Dunnivant  and 
Schwarzenbach  (1992)  reported  that  TNT  degradation  caused  by  natural  organic 
matter  (NOM)  was  increased  by  elevated  pH.  Previous  work  in  our  laboratory 
has  shown  that  degradation  rates  of  TNT  associated  with  base-induced  transfor¬ 
mation  technology  is  rapid,  with  complete  TNT  degradation  realized  in  <40  min 
at  room  temperature  (Felt,  Larson,  and  Hansen  2001a).  Kinetics  studies  of  TNT 
and  other  nitroaromatics  after  alkaline  hydrolysis  demonstrated  similar  results  in 
aqueous  solutions  and  in  highly  contaminated  soils  using  calcium  hydroxide 
(Emmrich  1999,  2001). 

Previous  kinetics  studies  concentrated  on  the  degradation  of  the  parent  com¬ 
pound,  TNT,  but  gathered  little  data  on  the  kinetics  of  the  over-all  OH-TNT 
reaction  mixture  as  a  whole.  There  is  also  a  need  to  determine  information  on 
the  transformation  products  of  the  OH-TNT  reaction.  A  common  analytical 
method  used  to  identify  individual  components  in  a  reaction  mixture  is  liquid 
chromatography-mass  spectroscopy  (LC-MS),  which  separates  the  components 
of  a  mixture  on  the  chromatography  column  and  identifies  them  by  comparing 
their  mass  spectrum  to  spectrum  of  known  standards  (Willard  et  al.  1988).  From 
a  previous  study  in  our  laboratory,  it  was  determined  that  roughly  50  percent  of 
the  final  products  of  the  OH-TNT  reaction  are  polymers,  with  molecular  weights 
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above  1,000  Daltons  (Felt,  Larson,  and  Hansen  2001b).  Polymers  cannot  be 
identified  using  LC-MS,  because  their  large  molecular  size  causes  them  to  foul 
the  chromatography  column.  An  alternative  method  that  could  yield  kinetic  data 
of  the  overall  reaction  and  theoretically  separate  the  reaction  components  is  a 
factor  analysis  of  UV/VIS  spectral  data. 

UV/VIS  spectral  analysis  is  based  on  the  Lambert-Beer  Law  that  relates  the 
absorbance  of  a  chemical  species  to  its  concentration, 

Ax  =  £x  bC  (1) 

where 

Ax  =  absorbance  of  a  compound  at  a  given  wavelength  (X) 

e  =  molar  absorptivity  of  the  compound  at  wavelength  X 

b  =  pathlength  of  the  sample  cell 

C  =  concentration  of  the  compound  (Willard  et  al.  1988) 

In  a  spectral  experiment,  a  real  experimental  data  matrix  of  absorbance  over  time 
[A\t)\  may  be  rationally  constructed  from  a  row  matrix  [£*]  of  molar  absorptivi- 
ties  and  a  column  matrix  [c(0]  such  that 

[A/)]  =  [ex][c(r)]  (2) 

The  pathlength,  b ,  of  the  sample  cell  is  usually  1  cm  and  it  remains  constant 
over  the  experiment,  so  it  is  dropped  from  Equation  2.  In  application,  there  will 
be  as  many  rows  in  [£*]  as  there  are  spectral  wavelengths  observed  and  as  many 
columns  as  there  are  reaction  components,  i  (Malinskowski  1989).  Similarly, 
there  will  be  in  [c(f)]  as  many  rows  as  there  are  reaction  components,  /,  and  as 
many  columns  as  there  are  times  observed  during  the  course  of  the  reaction. 

While  the  overall  dimensionality  for  the  system  (r  x  c)  is  known  from  the  design 
of  the  experiment,  the  number  of  reaction  components  is  hidden  in  the  matrix 
multiplication  where  an  [r  x  /]  matrix  is  multiplied  by  an  [/  x  c]  matrix  to  produce 
the  [r  x  c]  data  matrix.  The  minimum  number  of  abstract  terms  («)  required  to 
account  for  the  missing  reaction  component  information  can  be  determined 
without  a  priori  information  about  /  by  using  principal  component  analysis 
(PCA). 

The  real  data  matrix  and  its  abstract  row  and  column  matrix  constituents  are 
defined  as 

[D]  =  [/?][q  (3) 


The  raw  experimental  data  are  not  analyzed  directly  but  converted  to  a 
covariance  matrix  [Z]  by  multiplying  the  data  matrix  by  its  transpose 
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[Z]  =  [D]t[D] 


(4) 


The  diagonalization  of  this  matrix  by  the  eigenvalue-eigenvector  method 
produces 


[Z]  =  [C]T[X][C] 


(5) 


where  [X]  is  an  i  x  i  (diagonal)  matrix  of  eigenvalues  which  are  ordered  from 
most  to  least  important.  Only  the  first  and  larger  n  terms  of  [X]  account  for  these 
data,  with  the  remaining  i-n  terms  being  essentially  nil  and  attributable  to  noise. 
The  abstract  matrix  [C*]  is  the  diagonalization  matrix,  since 


[c*r=[c*r’ 

(6) 

[&W.cxV  =  [X] 

(7) 

Hence  [X\  is  essentially  an  nxn  matrix,  and  the  abstract  matrix  [C*]  contains  the 
eigenvectors  and  has  dimensions  nxc. 

The  abstract  row  matrix  [R+]  of  eigenvalues  can  then  be  calculated  from 

[Rl]  =  [D][C*]r  (8) 

which  has  rxn  dimensions.  These  data  matrix  can  then  be  reconstructed  from 
the  abstract  row  and  column  matrices  from  the  n  principal  components,  which  are 
derived  from  the  larger  terms  in  [X.]  using  PCA.  This  process  has  been  coded  into 
a  FORTRAN  program  “PCA.F90”  given  in  Appendix  A. 

PCA  produces  factors  which  are  implicitly  irrational,  orthogonal,  and 
account  for  [j D]  to  within  experimental  error.  The  chemical  details  of  the  system, 
a  rational  system,  are  embedded  in  [D]  and  cannot  be  revealed  by  a  simple 
inspection  of  the  eigenvalues  and  eigenvectors  of  PCA.  The  abstract  factors  can 
be  mathematically  transformed  into  rational  factors  by  using  target  testing  analy¬ 
sis  which  tests  suspected  factors  representing  a  row  of  the  row  matrix  or  column 
of  the  column  matrix  against  the  experimental  data  to  determine  if  it  is  a  real 
factor.  The  most  readily  testable  target  factors  are  initial  and  final  spectral  factors 
taken  directly  from  the  data.  Other  factors  may  suggest  themselves  on  the  basis 
of  the  kinetics  or  intermediates  observed  during  the  course  of  the  reaction. 

A  row  test  vector  R  "  (eigenvalue)  can  be  shown  using  a  least  squares  method 
to  be  a  real  factor  if  the  transformation  vector,  T'r,  implied  by 

T'r=[X]-l[Rt]R"  (9) 

produces  vector  7?'  which  closely  agrees  with  R"  on  application  as  follows: 
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R"~R'=[RX]  T'r 


(10) 


A  similar  target  approach  with  the  eigenvectors  in  [C*]  also  can  be  done.  In 
this  case,  a  single  column  test  vector  C"  (eigenvector)  is  a  real  factor  if  the  trans¬ 
formation  vector,  T'c,  implied  by 

rc  =  C"[Z][C*]W  (11) 


produces  vector  C'  which  closely  agrees  with  C"  on  application  as  follows: 

C" *  C' -  T'c[Cx]  (12) 

Individual  vectors  can  be  tested  in  either  the  row  or  column  categories  and 
discover  whether  they  are  real  factors  independent  of  any  undetermined  factors. 
By  doing  so,  and  developing  a  theoretical  or  experimental  model,  one  can  in 
principle  reconstruct  a  set  of  n  real  factors  that  accounts  for  the  data  to  within 
experimental  error. 

In  summary,  experimental  spectral  data  (absorbance  versus  time)  can  be  used 
to  create  a  data  matrix  of  molar  absorptivities  and  concentrations  containing 
information  about  each  reaction  component  in  the  reaction  mixture  over  time. 

The  row  matrix  in  the  kinetic  spectroscopic  application  contains  the  spectra  of 
the  n  components,  while  the  column  matrix  contains  the  time  relative  concentra¬ 
tion  data  for  the  same  n  components.  In  principle,  the  spectra  and  rate  data  for 
each  component  can  be  retrieved  using  PCA  and  target  testing  analysis  and  be 
used  to  develop  kinetic  models  of  the  reaction  using  a  least  squares  method 
without  the  use  of  separation  techniques  or  mass  spectrometry. 
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3  Materials  and  Methods 


Previous  experiments  had  studied  the  kinetics  of  the  degradation  of  TNT 
after  base  addition,  but  the  overall  reaction  kinetics  had  not  been  studied.  In  this 
study,  UV/VIS  spectral  analysis  of  the  OH-TNT  reaction  was  conducted  over 
time  and  at  different  temperatures.  An  FA  program  was  used  to  analyze  the 
spectral  data  and  determine  the  number  of  major  chemical  species  in  the  overall 
reaction  without  using  separation  technologies  (e.g.,  gas  chromatography  (GC), 
high-pressure  liquid  chromatography  (HPLC),  or  liquid  chromatography-mass 
spectroscopy  (LC-MS).  Test  spectral  vectors  were  developed  and  tested  against 
abstract  vectors.  Results  were  used  to  indicate  possible  kinetic  models. 


Materials 

Chemicals  and  glassware 

Chemicals  used  in  this  study  included  TNT  supplied  by  the  Rock  Island 
Arsenal,  Rock  Island,  IL,  reagent  grade  potassium  hydroxide  purchased  from 
Fisher  Scientific,  and  dionized  water.  Glassware  used  in  this  study  included 
beakers,  pastuer  pipettes,  and  volumetric  flasks. 


Instrumentation 

The  UV/VIS  spectrometer  was  a  Hewlett  Packard  8453  with  a  diode  array 
detector  (DAD)  with  a  1  -nanometer  (nm)  resolution  purchased  from  Agilent 
Technologies.  The  instrument  was  equipped  with  UV/VIS  HPChem  software  and 
a  jacketed  1.0-cm  quartz  sample  cell  that  was  thermostated  by  connection  to  a 
recirculating  water  bath  to  maintain  a  given  temperature. 


Methods 

Spectroscopy  experiment 

UV/VIS  spectral  analyses  of  the  OH-TNT  reaction  were  conducted  at  four 
different  temperatures  to  develop  information  concerning  the  reaction  mechanism 
and  the  individual  components  of  the  OH-TNT  reaction  as  they  change  over 
time. 
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Two  milliliters  (mL)  of  water  was  added  to  1  mL  TNT  solution  (102  ppm)  in 
a  25-mL  flask  and  the  solution  was  equilibrated  at  25  °C.  At  time  zero,  1  mL  1  N 
KOH  was  added  and  the  mixture  was  homogenized  by  swirling  briefly  to  yield  a 
reaction  mixture  of  2:1:1  (water:  100  ppm  TNT:  IN  KOH).  A  representative 
sample  (1  mL)  was  pipetted  into  a  jacketed  1 .0-cm  quartz  cell  and  analyzed  using 
UV/VIS  spectrometry  from  190  to  1,100  nm.  The  cell  was  thermostated  by  con¬ 
nection  to  a  recirculating  water  bath  maintained  at  25  °C.  These  data  were 
collected  in  a  darkened  room  to  prevent  TNT  degradation  as  a  result  of  photoly¬ 
sis.  A  diluted  (1 :3)  sample  of  the  initial  TNT  solution  was  analyzed  to  provide  a 
reference  spectrum. 

From  a  previous  kinetics  study,  it  was  known  that  the  OH-TNT  reaction  had 
a  very  fast  initial  decay  and  a  much  slower  secondary  decay  (Felt,  Larson,  and 
Hansen  2001a).  The  sample  spectral  runs  were  set  up  in  two  parts  to  take  this  into 
account.  The  spectrometer  was  programmed  to  yield  numerous  data  points  early 
in  the  reaction  to  measure  the  very  fast  initial  decay.  Allowance  for  a  longer 
sampling  cycle  during  the  second  part  of  the  test  is  necessary  because  the  second¬ 
ary  degradation  is  much  slower.  The  initial  settings  for  the  instrument  included  a 
10-  to  12-s  delay,  a  20-s  sampling  cycle,  and  a  total  run  time  of  1,200  s  (20  min). 
No  incremental  increase  in  the  sampling  cycle  was  programmed  for  this  part  of 
the  experiment.  Data  collection  settings  were  changed  after  the  first  protocol  to 
remove  the  delay  and  modify  the  sampling  cycle  to  300  s  with  an  increment  of 
25  percent  after  1,800  s,  with  an  additional  run  time  of  43,200  s  (12  h).  The 
instrument  was  started,  and  the  time  elapsed  between  the  two  runs  was  noted. 

The  experimental  protocol  described  above  was  repeated  using  incubation 
temperatures  of  20,  15,  and  12  °C. 


Factor  analysis 

FA  was  applied  to  the  experimental  data  using  the  method  described  by 
Malinskowski  (1989).  The  PCA  process  has  been  coded  into  a  FORTRAN  pro¬ 
gram  “PCA.F90”  and  is  given  in  Appendix  A.  The  program  used  to  test  vectors 
derived  from  PCA  is  given  in  Appendix  B. 

The  FA  method  selected  for  evaluation  is  a  multivariate  analysis  that  can  be 
applied  to  a  large  data  set  from  superposed,  multiple  sources  such  as  individual 
compounds  in  a  reaction  mixture  (Malinskowski  1989).  Chemical  transforma¬ 
tions  result  in  varying  quantities  of  reactants,  intermediates,  and  products  over 
time,  each  with  their  characteristic  spectral  signatures.  During  spectral  analysis,  a 
detector  measures  transmitted  light  intensity  (which  is  converted  to  the  derivative 
unit  “absorbance”)  over  a  range  of  wavelengths  and  at  time  intervals  during  the 
course  of  the  reaction.  Each  chemical  species  is  independently  associated  with  its 
own  absorbance  and  the  absorbance  of  the  reaction  mixture  that  is  indicated  by 
the  instrument  is  a  combination  of  these  signals.  Chemical  interference  is  possi¬ 
ble,  but  it  is  rare  and  was  assumed  to  be  inconsequential  in  this  study. 

In  the  FA  procedure  used  in  this  study,  the  experimental  data  are  represented 
by  a  spectral  data  matrix,  consisting  of  molar  absorptivities  and  concentrations  as 
a  function  of  time.  The  experimental  data  matrix  is  mathematically  converted  to 
an  abstract  matrix  that  is  analyzed  using  PCA  to  yield  the  minimum  number  of 
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factors  necessary  to  reconstruct  the  abstract  row  matrix  within  the  experimental 
error.  The  principal  components  are  irrational  factors  that  do  not  directly  repre¬ 
sent  the  actual  reaction  components  of  the  OH-TNT  reaction.  The  irrational 
factors  indicated  by  PCA  are  transformed  into  rational  (real)  components  using 
target  testing.  Test  spectral  vectors  for  four  factors,  including  the  reactant  (TNT), 
two  intermediates,  and  the  final  product,  were  developed  and  tested  against  the 
abstract  vectors.  There  was  good  agreement  between  the  test  vectors  and  the 
abstract  vectors,  indicating  the  test  vectors  were  real  factors,  corresponding  to 
real  components  of  the  reaction.  Two  kinetic  models  that  incorporated  the  test 
vectors  were  developed  and  refined  using  a  nonlinear  regression  method. 

Figure  1  is  a  flowchart  that  summarizes  the  FA  process. 
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Experimental  data 


Reconstruct  abstract  row  matrix 

i 

Construct  test  spectral  vectors  using  observations  of  experimental  data 

i 

Target  testing 

I 


Figure  1 .  Flowchart  for  factor  analysis  of  spectral  kinetic  data 
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4  Results 


Spectral  Data  Collected  at  25  °C 

Initial  TNT  spectrum 

UV/VIS  spectroscopy  charted  the  chemical  changes  that  occurred  during  the 
OH-TNT  reaction  from  the  initial  TNT  spectra  (Figure  2)  through  the  intermedi¬ 
ates  until  final  products  were  formed.  The  initial  TNT  spectrum  indicated  that 
TNT  absorbs  in  the  ultraviolet  range,  with  a  shoulder  at  254  nm.  The  spectrum  is 
flat  in  the  visible  range  (400  to  700  nm),  which  is  consistent  with  a  colorless 
fluid. 


Figure  2.  Spectrum  of  diluted  (1:3)  TNT  stock  solution 

Spectra  of  reaction  mixture  at  25  °C 

Figure  3  shows  the  spectra  generated  during  the  initial  run  at  25  °C.  The 
arrows  indicate  the  direction  of  the  spectral  changes  over  the  course  of  data 
collection.  The  first  spectrum  at  10  s  shows  a  defined  shoulder  at  240  nm.  The 
spectrum  reaches  a  minimum  at  360  nm,  rises  to  another  maximum  of  about 
0.8  AU  at  450  nm,  and  reaches  a  final  minimum  at  650  nm.  No  absorbance  is 
seen  from  650  to  1,100  nm.  The  second  and  subsequent  spectra  show  that  the 
240-nm  shoulder  has  been  shifted  to  260  nm.  The  minimum  has  shifted  to  340 
from  360  nm,  and  the  maximum  at  450  nm  has  shifted  slightly  to  the  left  and 
continues  to  rise.  This  maximum  reaches  its  peak  by  191  s,  before  it  slowly 
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Figure  3.  Spectra  of  OH-TNT  reaction  at  25  °C  -  first  20  min 

begins  to  decrease  over  the  remainder  of  the  reaction.  A  slight  shoulder 
developed  at  510  nm,  and  a  foot  developed  between  510  and  550  nm. 

Figure  4  shows  the  spectra  generated  during  the  latter  portion  of  the  experi¬ 
ment  at  25  °C.  The  most  obvious  feature  is  the  slow  decrease  in  the  spectral 
feature  at  450  nm.  The  spectra  have  an  approximate  isosbestic  point  at  340  nm. 
The  feature  at  330  nm  has  risen  and  remains  constant.  A  shoulder  is  still  obvious 
at  260  nm  and  slowly  decreases.  A  foot  is  still  visible  between  510  and  550  nm. 


Spectrum  of  Final  Products 

The  spectrum  of  the  final  products  was  collected  after  48  h  at  which  time  the 
reaction  was  judged  to  be  complete  and  is  illustrated  in  Figure  5.  A  broad  feature 
at  330  nm  and  a  shoulder  at  270  nm  are  evident.  The  spectrum  has  a  low 
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Figure  4.  Spectra  of  OH-TNT  reaction  at  25  °C  - 1 2  h 


Figure  5.  Spectrum  of  final  products  of  OH-TNT  reaction 
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absorption  in  the  visible  range  decreasing  for  400  to  600  nm.  This  result  is  con¬ 
sistent  with  the  light  straw  yellow  color  of  the  final  reaction  mixture. 


Spectral  data  collected  at  20, 15,  and  12  °C 

In  order  to  collect  more  kinetic  data  about  the  OH-TNT  reaction,  the  experi¬ 
ment  described  above  was  repeated  at  different  temperatures  (20,  15,  and  12  °C). 
These  spectra  showed  the  rise  of  a  major  spectral  feature  with  a  maximum  of 
0.8  AU  at  450  nm  at  all  temperatures,  but  the  features  appear  after  longer  times  at 
cooler  temperatures.  These  data  correspond  to  the  data  recorded  at  25  °C,  except 
that  the  maximum  was  reached  more  slowly  at  the  lower  temperature  (191  s  at 
25  °C,  230  s  at  20  °C,  and  419  s  at  15  °C). 

Factor  analysis  results 

PCA  indicated  that  six  principal  components  explained  the  spectra  to  within 
experimental  error,  with  four  factors  explaining  the  majority  of  the  variance.  Test 
spectral  vectors  for  four  components  were  developed,  including  TNT,  two  inter¬ 
mediates,  and  the  final  product,  and  were  tested  against  abstract  vectors.  There 
was  good  agreement  between  the  test  vectors  and  the  abstract  vectors,  indicating 
the  test  vectors  were  real  factors,  corresponding  to  real  components  of  the 
reaction. 
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Spectroscopy 

Temperature  dependence  of  OH-TNT  reaction 

A  previous  study  using  HPLC  indicated  that  the  second  phase  of  the  OH- 
TNT  reaction,  the  slower  secondary  degradation,  is  temperature-dependent  (Felt, 
Larson,  and  Hanson  2001a).  The  rate  of  this  phase  of  the  OH-TNT  reaction 
decreased  as  the  temperature  decreased.  In  this  study,  it  appears  that  the  OH-TNT 
reaction  proceeded  by  the  same  mechanism  at  all  temperatures  that  were  moni¬ 
tored  and  the  differences  in  the  spectra  were  a  result  of  the  temperature 
dependence  of  the  OH-TNT  reaction. 

The  450-nm  feature  reached  a  maximum  of  about  0.8  AU,  regardless  of  the 
incubation  temperature  and  the  run  time  that  the  feature  occurred.  This  may  indi¬ 
cate  that  the  450-nm  feature  is  a  pure  spectral  feature  of  the  second  intermediate 
and  that  this  component  has  a  large  concentration  and  a  relatively  small  molar 
absorptivity. 


Primary  and  secondary  intermediates 

There  was  spectral  and  visual  evidence  that  indicated  a  primary  intermediate 
was  formed  and  quickly  replaced  by  a  secondary  intermediate.  In  the  current 
study,  the  spectral  feature  at  450  nm  on  the  first  spectrum  (Al)  at  25  °C  was 
visually  different  in  shape  than  those  of  successive  spectra  at  that  incubation 
temperature.  The  reaction  mixture  changed  from  a  colorless  solution  to  a  pink 
shade  after  addition  of  the  base,  and  then  quickly  changed  to  a  reddish  orange, 
indicating  at  least  two  separate  reaction  intermediates.  The  OH-TNT  reaction 
mixture  was  a  light  yellow  color  at  the  end  of  48  h,  indicating  the  final  product  is 
another  chemical  species. 

In  order  to  determine  if  a  primary  intermediate  had  been  formed,  Al  was 
scaled  and  compared  to  a  successive  spectrum  at  25  °C.  If  Al  contained  infor¬ 
mation  associated  with  a  different  component  than  the  successive  spectra  at 
25  °C,  Al  would  have  a  different  shape  than  that  of  the  other  spectra.  If  the 
components  of  Al  and  A10  were  the  same,  the  two  peaks  should  be  different 
from  each  other  only  in  size,  not  shape.  Al  was  scaled  with  the  spectra  that 
yielded  the  450-nm  maximum  (A  10,  time  =  191  s)  and  is  illustrated  as  the  thin 
solid  line  in  the  middle  of  Figure  6.  The  scale  was  the  absorbance  of  A10  at 
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Wavelength  (nm) 


Figure  6.  Spectral  features  of  the  first  intermediate  of  OH-TNT  reaction  at 
25  °C 


450  nm  divided  by  the  absorbance  of  A1  at  450  nm.  A  plot  of  the  difference  of 
scaled  A1  less  A10  in  absorbance  units  versus  wavelength  would  be  a  straight 
line  through  zero  if  the  reaction  components  that  make  up  A1  were  the  same  as 
those  that  make  up  A 10.  When  the  difference  was  plotted  (the  bold  line  in  Fig¬ 
ure  6),  the  line  did  not  go  through  zero  but  increased  to  a  maximum  at  550  nm. 
This  could  indicate  that  a  reaction  component  appears  at  11  s  (Al),  which  is 
distinct  from  the  species  producing  the  large  450-nm  peak  that  reaches  its  maxi¬ 
mum  at  191  s  (A  10).  The  first  derivative  of  the  difference  also  shows  inflections 
at  450  and  490  nm,  as  shown  by  the  dashed  line  in  Figure  6.  Spectrum  Al  is  also 
different  than  the  initial  TNT  spectra,  which  would  indicate  an  intermediate  had 
formed  in  this  early  reaction  time.  This  intermediate  seemed  to  be  very  reactive, 
as  the  absorbance  of  500  nm  at  32  s  (A2)  had  already  decreased  relative  to  the 
first  spectra.  This  indicates  that  a  second  intermediate  had  already  started  to  form 
by  32  s  into  the  reaction,  which  appeared  to  be  producing  the  significant  absorp¬ 
tion  at  450  nm. 

These  conclusions  are  consistent  with  the  color  changes  noted  for  the  reac¬ 
tion  mixtures  and  results  from  a  previous  kinetic  study  of  TNT  degradation  after 
introduction  of  base  using  HPLC  (Felt,  Larson,  and  Hansen  2001a).  The  results 
of  this  study  indicated  a  two-phase  reaction,  a  very  fast  initial  rate,  followed  by  a 
much  slower  degradation  rate,  which  is  consistent  with  the  formation  and  degra¬ 
dation  of  two  intermediates.  Chromatograms  of  the  base-challenged  TNT  solu¬ 
tions  showed  a  feature  that  had  formed  quickly  at  a  3.5-min  retention  time  that 
was  indicative  of  an  intermediate.  The  peak  height  of  this  intermediate  dimin¬ 
ished  over  time,  indicating  follow-on  degradation.  The  short  retention  time  on  the 
reverse  phase  C-l  8  column  indicated  a  compound  that  was  more  polar  than  TNT 
(retention  time  17.5  min).  The  majority  (98  percent)  of  the  final  reaction  products 
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of  the  OH-TNT  reaction  are  also  polar  compounds,  as  indicated  by  a  study  that 
separated  the  reaction  products  by  molecular  weight  and  solubility  (Felt,  Larson, 
and  Hansen  2001b). 


Factor  Analysis  of  Spectral  Data 

Chemical  transformations  often  result  in  varying  quantities  of  starting  mate¬ 
rial,  intermediates,  and  products  over  time,  each  with  their  characteristic  spectral 
signatures.  FA  has  particular  applicability  in  these  chemical  systems  in  which 
detector  output  records  events  containing  information  from  superposed,  multiple 
components.  Absorbance  was  assumed  to  be  independently  associated  with  each 
chemical  species.  Chemical  interference  was  possible,  but  it  was  assumed  to  be 
inconsequential. 

Absorbance  data  for  multicomponent  mixtures  can  be  approximated  using 
the  Lambert-Beer  Law,  where  each  data  point  can  be  represented  by  the  equation, 
Aik=  X  G  ijCjk-  s  ij  is  the  molar  absorptivity  per  unit  path  length  of  component  j  at 
wavelength  i,  and  cjk  is  the  concentration  of  component  j  in  the  kth  mixture. 
Spectral  data  can  theoretically  be  explained  by  a  hypothetical  experimental  data 
matrix  containing  row  matrix  [e^]  of  molar  absorptivities  (spectra)  and  column 
matrix  [c(t)  ]  of  time-dependent  concentrations  (Malinskowski  1989).  FA  identi¬ 
fies  the  number  of  principal  reaction  components  in  the  reaction  mixture  and 
produces  irrational,  orthogonal  factors  that  account  for  this  experimental  matrix 
to  within  experimental  error.  The  abstract  matrix  can  be  converted  into  real 
factors  (chemically  relevant  data)  by  either  using  an  appropriate  transformation 
matrix  or  by  using  target  testing,  which  is  the  method  used  in  this  study.  In  order 
to  test  the  evolving  model,  suspected  factors  must  be  identified  that  can  be  used 
as  starting  points  for  the  model.  Factors  are  suggested  on  the  basis  of  spectral  and 
other  experimental  data  and  tested  against  the  spectral  data  to  determine  if  they 
are  real  factors. 

PCA  indicated  six  factors  were  required  to  reduce  the  reproduced  data  matrix 
to  within  the  estimated  measurement  errors.  Using  four  principal  components 
leaves  a  residual  data  matrix  with  all  elements  less  than  0.010  AU;  the  estimated 
measurement  errors  were  less  than  this,  at  0.002  AU.  The  two  additional  factors 
were  very  small  contributors  to  the  spectral  variance,  were  without  spectral 
features  of  note,  and  were  therefore  ignored.  Table  1  shows  the  results  of  the 
PCA  of  the  25  °C  spectral  data.  Eigenvalues  are  listed  in  order  of  the  largest  to 
the  smallest,  and  the  associated  eigenvector  for  each  time  is  given.  The  size  of 
the  eigenvalue  signifies  the  amount  of  the  variance  in  these  data  that  is  explained 
by  that  eigenvalue.  The  eigenvectors  associated  with  the  eigenvalues  mimic  con¬ 
centration  versus  time  for  the  abstract  components.  The  first  eigenvalue,  517, 
indicates  that  this  eigenvalue  explains  most  of  the  spectral  data.  The  next  three 
eigenvalues  are  significant  (22.47,  1.02,  and  0.27)  and  represent  the  other  three 
principal  components  of  the  reaction  as  identified  by  PCA. 
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Table  1 

Spectral  Vectors  for  TNT-Hydroxide  Reaction  at  25  °C 
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0.14 

0.12 

0.09 

0.02 

j  -0.06  -0.13  -0.17 

-0,16  -0.11 

-0.03 

-0.01 

-0.01 

-0.02 

-0.01 

-0.01  0.10  0.22 

0.31 

REAL  ERROR-  0.00  EMBEDDED  ERROR-  0.00 

EIGENVALUE  ft  8 

0.001 

-0.12  0.57  -0.41 

-0.32  -0.13 

-0.01 

0.05 

0.08 

0.12 

0.10 

0.10  0.06  0.04  . 

0.00  .  0.02 

-0.07 

-0.04 

-0.07 

-0.02 

-0.04 

-0.04  0.05  0.04 

0.89  0.02 

0.04 

0.06 

0.04 

0.03 

0.05 

-0.02  -0.09  -0.12 

-0.14  -0.17 

-0.10 

-0.08 

-0.07 

-0.07 

-0.07 

-0.04  0.10  0.21 

0.37 

REAL  ERROR-  0.00  EMBEDDED  ERROR-  0.00 

EIGENVALUE  ft  9 

0.000 

0.00  -0.08  0.26 

0.07  -0.07 

-0.08 

-0.13 

-0.09 

-0.14 

-0.12 

-0.12  -0.02  0,05 

0.06  0.08 

0.24 

0.19 

0.19 

0.12 

0.18 

0.07  -0,14  -0.18 

-0.26  -0.17 

-0.15 

-0.06 

0.04 

0.23 

0.23 

0.14  -0.03  -0.14 

-0.16  -0.07 

-0.09 

-0.09 

-0.12 

-0.12 

-0.06 

-0.10  -0.03  0.25 

0.41 

REAL  ERROR-  0.00  EMBEDDED  ERROR-  0.00 

EIGENVALUE  ft  10 

0.000 

1  -0.02  0.12  0.01 

-0.16  -0.13 

-0.02 

0.04 

.0,07 

-0.01 

-0.02 

|  -0.02  0.10  0.00 

0.00  -0.03 

0.26 

0.09 

0.14 

-0.05 

0.03 

-0.09  -0,30  -0.34 

-0.20  0.25 

0.25 

0.23 

0.21 

0.07 

-0.25 

-0.22  -0.20  -0.19 

-0.07  0.18 

0.18 

0.15 

0.13 

0.04 

0.11 

-0.02  -0,14  -0.07 

-0.11 

REAL  ERROR-  0,00  EMBEDDED  ERROR-  0.00 
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The  eigenvalues  and  eigenvectors  were  used  to  reconstruct  the  abstract  row 
matrix.  Row  matrices  (Table  2)  mimic  the  spectral  data,  are  part  of  an  abstract 
matrix,  but  may  have  little  rational  significance.  The  four  factors  listed  in  Table  2 
represent  the  four  principal  components  identified  by  PCA,  as  described  above. 
The  listed  values  illustrate  trends  of  the  abstract  matrix  compared  to  measurable 
spectral  values. 


Table  2 

Reconstructed  Abstract  Row  Matrix  for  Four  Important  Factors 
for  OH-TNT  Reaction  at  25  °C 

Wavelength  (nm) 

1 

2 

3 

4 

220. 

11.40 

1.81 

0.02 

0.06 

230. 

7.44 

1.11 

0.25 

0.05 

240. 

6.80 

0.72 

0.25 

-0.01 

250. 

6.10 

0.57 

0.21 

-0.04 

260. 

5.22 

0.40 

0.17 

-0.06 

270. 

4.59 

0.23 

0.09 

-0.09 

280. 

4.17 

0.10 

-0.08 

-0.12 

290. 

3.49 

0.09 

-0.19 

-0.07 

300. 

2.93 

0.11 

-0.20 

-0.02 

310. 

2.55 

0.16 

-0.23 

0.01 

320. 

2.27 

0.25 

-0.27 

0.02 

330. 

2.10 

0.34 

-0.30 

0.03 

340. 

2.00 

0.39 

-0.32 

0.04 

350. 

1.95 

0.36 

-0.33 

0.04 

360. 

1.90 

0.24 

-0.31 

0.04 

370. 

1.97 

0.06 

-0.28 

0.02 

380. 

2.08 

-0.18 

-0.23 

0.00 

390. 

2.28 

-0.42 

-0.18 

-0.01 

400. 

2.55 

-0.67 

-0.14 

-0.03 

410. 

2.86 

-0.89 

-0.09 

-0.05 

420. 

3.15 

-1.08 

-0.05 

-0.07 

430. 

3.35 

-1.21 

-0.01 

-0.09 

440. 

3.55 

-1.34 

0.02 

-0.10 

450. 

3.58 

-1.38 

0.05 

-0.10 

460. 

3.51 

-1.38 

0.06 

-0.07 

470. 

3.31 

-1.32 

0.06 

-0.02 

480. 

3.00 

-1.20 

0.05 

0.05 

490. 

2.65 

-1.04 

0.04 

0.12 

500. 

2.30 

-0.87 

0.04 

0.15 

510. 

2.01 

-0.73 

0.03 

0.15 

520. 

1.75 

-0.62 

0.03 

0.14 

530. 

1.48 

-0.51 

0.03 

0.13 

540. 

1.24 

-0.42 

0.04 

0.14 

550. 

1.04 

-0.34 

0.04 

0.14 

560. 

0.88 

-0.28 

0.03 

0.13 

570. 

0.75 

-0.22 

0.02 

0.11 

580. 

0.64 

-0.17 

0.01 

0.09 

590. 

0.54 

-0.14 

0.01 

0.06 

600. 

0.47 

-0.11 

0.00 

0.05 

610. 

0.41 

-0.08 

-0.01 

0.03 

620. 

0.38 

-0.06 

-0.01 

0.02 

630. 

0.34 

-0.05 

-0.01 

0.02 

640. 

0.31 

-0.03 

-0.01 

0.01 

650. 

0.29 

-0.03 

-0.01 

0.01 

660. 

0.27 

-0.02 

-0.01 

0.01 

670. 

0.26 

-0.02 

-0.01 

0.01 

680. 

0.25 

-0.01 

-0.01 

0.01 

690. 

0.25 

-0.01 

-0.01 

0.01 

700. 

0.24 

-0.01 

-0.01 

0.00 
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The  experimental  data  were  then  examined  for  information  that  would  sug¬ 
gest  possible  factors.  A  factor  for  one  component,  the  reactant,  was  taken  directly 
from  the  initial  TNT  spectrum.  The  first  few  spectra  of  the  OH-TNT  reaction 
revealed  a  possible  spectrum  of  a  material  with  a  broad,  possibly  divided  absorp¬ 
tion  centered  about  510  nm  (Figure  3).  This  was  rapidly  overcome  by  the 
appearance  of  a  species  with  a  strong  absorption  at  450  nm  that  built  to  a  maxi¬ 
mum  intensity  at  191  s,  then  slowly  declined  over  the  next  several  hours.  Late  in 
the  reaction,  a  chemical  species  with  a  significant  absorption  feature  at  330  nm 
was  noted  and  comprised  the  final  spectra. 

Test  spectral  vectors  were  constructed  for  the  second  intermediate,  the  final 
product,  initial  TNT,  and  the  first  intermediate  using  these  observations  of  the 
experimental  data.  The  second  intermediate  was  modeled  by  the  spectrum  for 
which  the  absorbance  at  450  nm  reaches  a  maximum  (191  s  at  25  °C)  and  corre¬ 
sponds  to  spectral  vector  1  on  Table  2.  The  values  start  off  very  high,  drop  to  a 
minimum  at  360  nm,  rise  to  a  maximum  at  450  nm,  and  then  continue  to  fall  to 
the  end,  which  corresponds  to  the  spectral  data.  The  final  product  was  modeled 
by  the  last  spectrum  (at  161,400  s  at  25  °C)  and  relates  to  spectral  vector  2.  The 
values  of  this  row  start  very  high,  drop  slightly,  come  back  to  a  maximum  at 
340  nm,  then  quickly  decrease  before  reaching  400  nm.  This  is  similar  to  the 
trend  seen  at  the  24-h  spectra  at  25  °C.  An  initial  TNT  spectrum  was  used  for  the 
third  spectral  vector.  This  row  begins  with  a  maximum  at  240  nm,  drops  quickly 
to  a  minimum  at  350  nm,  and  is  essentially  constant  at  higher  wavelengths.  Most 
of  the  features  of  TNT  are  in  the  UV  range,  with  a  maximum  at  240  nm,  and  the 
spectrum  is  constant  in  the  visible  range  (400  to  700  nm).  The  fourth  vector  was 
modeled  by  altering  the  TNT  spectrum  with  a  broad  double  peak  with  absorbance 
of  0.1  AU  centered  at  520  nm,  which  corresponds  to  spectral  vector  4.  This  row 
exhibits  a  maximum  at  340  nm  and  a  divided  peak  at  510  and  550  nm,  which  is 
consistent  with  the  spectral  data.  Overall,  there  is  good  agreement  between  the 
trends  of  the  factors  in  the  reconstructed  row  matrix  and  trends  in  the  observed 
spectra  of  the  reaction  over  time. 

The  model  that  emerges  is  a  four-component  reaction  with  the  parenthetical 
distinctive  features: 

TNT  ->  Intermediate  1  (530  nm)  -^Intermediate  2  (450  nm) 

Final  Products  (340  nm) 

Fast  Very  fast  Slow 


Test  Vectors  for  Four  Components 

Malinskowski  has  shown  that,  given  an  FA  solution,  the  abstract  row  vectors 
can  be  transformed  into  a  matrix  of  real  (chemically  meaningful)  row  vectors,  if 
the  transformation  matrix  can  be  established  and  if  it  describes  the  details  of  the 
transformation  of  the  row  vectors.  In  order  to  check  a  theoretical  model,  spectral 
vectors  must  be  identified,  as  described  above,  and  used  as  starting  points.  Table 
3  shows  the  abstract  row,  test,  and  confirmation  spectral  vectors  for  the  reaction 
at  25  °C  for  220-  to  700-nm  wavelengths.  The  vectors  marked  “rows”  were 
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Table  3 

Abstract  Row,  Test,  and  Confirmation  Spectral  Vectors  for  the  25  °C  Reaction 

1st  Spectral  Vector  SSI>=  0.00095 

ROW  1 1 .40  7.444  6.804  6.096  5.21 7  4.595  4.1 73  3.495  2.930  2.547  2.273  2.099  2.003  1 .946  1 .897  1 .968  2.079  2.279 

2.554  2.860  3.150  3.348  3.547  3.576  3.505  3.3  13  3.005  2.652  2.300  2.006  1.749  1.483  1.237  1.040  0.884  0.752  0.640  0.537 
0.466  0.4  13  0.377  0.340  0.306  0.287  0.273  0.262  0.255  0.247  0.243 

Trial  1 .632  1 .080  1 .033  0.937  0.81 1  0.729  0.673  0.565  0.471  0.400  0.341  0.297  0.277  0.274  0.286  0.327  0.382  0.453 

0.535  0.6  17  0.689  0.734  0.782  0.790  0.780  0.748  0.690  0.6  19  0.539  0.467  0.402  0.339  0.282  0.237  0.200  0.169  0.141  0.117 
0.099  0.084  0.074  0.065  0.057  0.052  0.048  0.046  0.044  0.042  0.041 

Confirmation  1.639  1.082  1.028  0.931  0.806  0.725  0.671  0.564  0.471  0.402  0.344  0.301  0.2790.2740.283  0.321  0.375  0.444 
0.526  0.609  0.684  0.735  0.788  0.799  0.789  0.753  0.689  0.6  0.532  0.461  0.400  0.339  0.284  0.240  0.203  0.170  0.142  0.117  0.098 
0.084  0.074  0.066  0.058  0.053  0.0500.04  0.045  0.044  0.043 

2nd  Spectral  Vector  SSE=  0.00274 

ROW  1.809  1.106  0.724  0,565  0.404  0.228  0.103  0.089  0.1100.1570.247  0.344  0.386  0.357  0.244  0.060  -.177  -.423 

-.667  -.888  -1.08  -1.21  -1. 34-1 .38-1 .38-1. 32-1.20-1 .04-.869-.726  -.615  -.512  -.420  -.343  -.279  -.222  -.175  -.137  -.105  -.079  -.058 
-.046-034-.028-.021-.017-.014-.012  -.01 1 

Trial  1 .829  1 .1 39  0.960  0.839  0.698  0.585  0.521  0.457  0.405  0.379  0.374  0.381  0.382  0.370  0337  0.302  0.2.59  0.226 

0.202  0.185  0.172  0.158  0.147  0.137  0.128  0.1190.113  0.106  0.099  0.093  0.083  0.073  0.063  0.056  0.051  0.048  0.044  0.039 
0.037  0.036  0.035  0.033  0.031  0.030  0.029  0.029  0.028  0.028  0.028 

Confirmation  1.807  1.130  0.964  0.847  0.707  0.601  0.545  0.476  0.416  0.384  0.376  0.380  0.381  0.369  0.335  0.299  0.253  0.217 
0.191  0.175  0.163  0.153  0.143  (7134  0.126  0.1180.1120.1070.103  0.0980.0900.078-.--67  0.058  0.053  0.049  0.045  0.040  0.038 
0.037  0.038  0.035  0.033  0.032  0.032  0.03  1  0.031  0.030  0.030 

3rd  Spectral  Vector  SSE=  0.01702 

ROW  0.018  0.249  0.250  0.205  0.167  0.092  -.083  -.186  -.203  -.225-.268-.302  -.324  -.329-.307  -.277  -.225  -.178  -.135 

-.095  -.054  -.0  14  0.022  0.045  0.056  0.055  0.048  0.042  0.036  0.030  0.027  0.032  0.039  0.039  0.034  0.023  0.014  0.006  -.002  -.007 
-.012  -.012  -.013  -.013  -.013  -.013  -.014  -.013  -.014 

Trial  1.759  1.332  1.166  1.004  0.821  0.623  0.374  0.214  0.155  0.121  0.092  0.083  0.074  0.060  0.044  0.030  0.018  0.010 

0.005  0.003  0.003  0.003  0.003  0.003  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002  0.002 
0.001  0.002  0.002  0.002  0.001  0.002  0.001  0.001  0.001  0.001  0.001 

Confirmation  1.763  1.315  1.141  0.987  0.817  0.644  0.427  0.262  0.187  0.140  0.102  0.085  0.069  0.049  0.028  0.007  -.009  -.022 
-.029  -.028  -.018  -.002  0.010  0.019  0.017  0.007  -.004  -.010  -.006  0.000  0.006  0.013  0.019  0.021  0.021  0.018  0.015  0.013  0.011 
0.011  0.011  0.011  0.011  0.011  0.011  0.011  0.011  0.011  0.011 

4th  Spectral  Vector  SSE  0.01078 

ROW  0.061  0.046  -.010  -.038  -.059  -.092  -.121  -.069  -.016  0.008  0.018  0.027  0.037  0.042  0.035  0.020  0.003  -.013  -.031 

-.050  -.070  -.088  -.101  -.098  -.074  -024  0.048  0.116  0.153  0.154  0.137  0.133  0.139  0.140  0.129  0.109  0.085  0.063  0.046  0.032 
0.023  0.017  0.012  0.010  0.008  0.007  0.006  0.005  0.005 

Trial  1.759  1.332  1.166  1.004  0.821  0.623  0.374  0.214  0.155  0.121  0.092  0.083  0.074  0.060  0.044  0.030  0.018  0.010 

0.005  0.003  0.003  0.003  0.003  0.006  0.010  0.025  0.040  0.080  0.100  0.100  0.095  0.098  0.100  0.100  0.080  0.060  0.040  0.020 
0.01 6  0.01 5  0.01 3  0.01 2  0.01 1  0.01 0  0.008  0.006  0.005  0.004  0.004 

Confirmation  1.779  1.337  1.147  0.981  0.803  0.616  0.383  0.236  0.180  0.140  0.103  0.085  0.071  0.054  0.033  0.011  -.004  -.015  - 
.022  -.021  -.013  0.000  0.011  0.024  0.032  0.041  0.056  0.075  0.089  0.091  0.086  0.088  0.093  0.093  0.086  0.073  0.058  0.045  0.034 
0.027  0.023  0.020  0.017  0.016  0.015  0.015  0.015  0.014  0.014 


the  same  rows  from  the  matrix  identified  in  Table  2.  The  vectors  marked  “Trial” 
are  the  test  spectral  vectors  constructed  previously  based  on  experimental  obser¬ 
vations.  The  rows  labeled  “Confirmation”  are  the  vectors  that  resulted  after  trans¬ 
formation  of  the  abstract  row  vectors  with  the  transformation  matrix.  Sum  of 
squares  error  (SSE)  in  this  method  is  defined  as: 

SSE  =  Z  (A1  (t)  ca,o  -  A\f)  obs)  (13) 

The  four  test  vectors  appear  to  agree  with  the  transformed  vectors  as  judged 
by  the  small  errors  between  them.  This  suggests  that  the  four  test  vectors  are  real 
factors  (spectral  components)  in  the  reaction. 
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Modeled  Spectra 

The  spectra  of  the  reaction  components  were  then  modeled  and  the  concen¬ 
trations  of  the  reaction  components  were  determined  over  time  using  the  test 
spectral  vectors  and  the  abstract  matrix.  Absorbance  data  of  base  hydrolysis  of 
TNT  were  interpreted  as  consistent  with  a  reactant,  two  intermediates,  and  a 
product  and  suggests  that  the  reaction  may  be  crudely  modeled  by  a  sequence  of 
first-order  reactions.  A  model  was  assembled  for  the  sequential  first-order 
reaction  scheme 

a  •¥  b  c  d 

The  composite  absorbance  for  this  reaction  scheme  at  any  time  and  wave¬ 
length,  Ay,  is  the  sum  of  the  absorbances  of  all  the  species  present  and  follows 
from  the  Lambert-Beer  Law, 

Ay  =  [a]ga  +  [b]gb  +  [c]ge  +  [d\gd  =  YL  gufikj  ( 1 4) 

where  the  g’s  are  the  molar  absorptivity  vectors  taken  over  the  /'  wavelengths, 
that  is,  they  are  the  spectra  of  the  absorbing  species,  and  the  concentrations  of  the 
k  species  are  followed  over  the  j  times.  Substituting  for  g’s  from  the  abstract  row 
vectors  Vk,  we  have 

A  =  [tf](TnLJ  +^12^2  ^13^3  +  ^14^4) 

+  [^]('C21^+  ^22^2  4*  ^23^3  +  X24V4) 

+  [c](  x3ivt  +  xnv2  +  x33V3  +  x34V4) 

+  MM  +  %2V2  +  x43V3  +  x43V4) 

and  collecting  terms  in  V,  we  have 

A  —  [(q)Tjj  +  (6)x21  +  (c)Tji  +  (cf)x4 j]  V{ 

+  [(a)x2l  +  (6) t22  +  (c)x32  +  (d)x42 ]  V2 
+  [(a)x31  +  (b)x32  +  (c)x33  +  (d)x34 ]  V3  ' 

+  [(fl)x41  +  (b)x42  +  (c)x43  +  (d)x44  ]  V4 

Defining  £  as  the  coefficients  of  V,  defined  over  the  1  components  as 


kkl  ~  In  Cl 


the  composite  spectrum  at  any  time  is 


Ai=Vl$u+V2Z>2l  +  ...  =  Vlk$k] 


(17) 


(18) 
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and  the  whole  absorbance  versus  time  data  set  becomes 


Au  —  ZZZ  VikTu-Cij  (19) 

It  is  seen  that  the  number  of  the  square  transformation  matrix  elements  (x’s) 
is  1  x  k,  where  1  =  k  and  the  number  of  important  factors  accounting  for  these 
data. 

In  order  to  derive  the  concentration  of  each  reaction  component  over  time, 
knowledge  of  the  rate  constants  ( ki ...  kk.\)  for  the  reaction  is  required.  For 
example,  the  concentration  of  the  k'h  component,  ck,  is 

ck(t)  =  CieV  +  C2ek2'...  +  Cke-kk‘  (20) 

where  C,  =  (kh  k2...k  k-i )  /  (k2-k})(k3-ki)...(kk-k]) 


and 


where  C2  =  (k\,  k2...k  k-\ )  /  (k\-k2)(k3-k2)...(kk-k2),  etc.  (Friedlander,  Kennedy,  and 
Miller  1964). 

Estimates  of  the  rate  constants  were  approached  by  ignoring  the  fast  initial 
reaction  (and  ignoring  the  first  few  data)  and  examining  the  behavior  of  the 
reaction  component  producing  the  450-nm  absorption.  Since  the  model  is  non¬ 
linear  in  all  of  the  variables,  the  variables  were  refined  against  the  absorbance 
versus  time  data  array  using  a  nonlinear  least  squares  procedure. 

Spectral  components  were  separated  by  a  nonlinear  least  squares  method. 
The  chi-squared  definition,  %2,  was  used  as  an  appropriate  gauge  of  the  agree¬ 
ment  between  a  developing  model  and  the  entire  spectra  when  fitting  the  whole 
spectra  measured  over  time.  The  best  model  for  k  components  would  have  a  set 
of  variables  that  produces  a  minimum  value  for  yj.  The  spectra  are  given  as  a 
function  of  the  PCA  row  vectors  and  a  transformation  matrix,  T 

(spectra)  uk  =u(Tmn)  (row)ik 
m-\...k\n  =  \...k 


and  the  concentrations  are  given  by  a  kinetic  model  tied  to  the  rate  constants. 


(concentrations)^-  =  u  (k\, ...  k  k.\  )kj 


(22) 


There  are  1[=  k2  +  (A>1)]  variables  for  k  components.  These  are  the  normal 
equations. 


“2  ^  [Wjy(^obsy  “  -^calc  ij  )][^4calc  ] ] 


(23) 
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In  this  case,  drflda \  is  a  vector  with  /  elements. 

In  order  to  obtain  the  best  fit,  the  variables  are  considered  parameters  to  be 
adjusted  to  minimize  %2-  The  derivative  of  is  taken  with  respect  to  each  vari¬ 
able;  i.e.,  the  gradient  of  %2  and  set  to  zero.  Taking  the  partial  derivatives  of  each 
of  the  observational  equations  with  respect  to  each  variable,  and  substituting  into 
the  normal  equations,  the  normal  equations  can  be  rearranged  to  segregate  the 
changes  in  the  variables,  which  minimize  %2.  These  equations  are  solved  by  a 
Gaussian  (linear)  elimination  method,  whereupon  the  shifts  are  applied  to  the 
variables,  and  the  process  is  repeated  until  convergence.  Because  the  system  is 
nonlinear,  the  shifts  will  often  result  in  sudden,  large  increases  in  rather  than  a 
smooth  procession  to  the  global  minimum.  Methods  to  handle  this  require  a 
combination  of  the  steepest  descent  method  (linear  least  squares)  with  dampening 
of  the  shifts  as  the  problem  moves  toward  the  minimum.  The  beginning  model 
must  resemble  the  final  solution,  since  the  polydimensional  hypersurface  is 
potentially  filled  with  many  false  minima  to  which  the  calculation  can  descend. 

After  the  variables  were  refined,  the  spectra  of  the  four  reaction  components 
were  predicted  by  the  sequential  first  order  kinetic  model  and  were  plotted  as 
illustrated  in  Figure  7.  The  four  components  in  this  figure  are:  TNT,  first  inter¬ 
mediate,  second  intermediate,  and  the  final  product.  TNT  spectra  plotted  in 
Figure  7  mirror  the  observed  spectra  in  Figure  2.  The  spectrum  shows  a  maxi¬ 
mum  at  250  nm,  decreases  quickly,  and  remains  featureless  in  the  visible  range. 


Wavelength  (nm) 


Figure  7.  Spectral  features  of  four  components  at  25  °C  using  sequential  first- 
order  kinetic  model 
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The  first  intermediate  shows  features  at  330  and  500  nm,  which  are  similar  to  the 
features  seen  growing  into  the  observed  spectra.  The  model  indicates  that  the  first 
intermediate  is  a  major  component  of  the  450-nm  feature  because  its  spectrum  is 
strong  at  that  wavelength.  Spectrum  of  the  second  intermediate  also  exhibits  a 
dominant  450-nm  feature  that  rises  quickly  and  decreases  slowly  over  the 
remainder  of  the  reaction,  comparable  to  the  observed  spectra.  The  spectrum  of 
final  products  calculated  by  this  model  exhibits  a  broad  feature  at  330  nm,  similar 
to  the  observed  spectra  of  the  final  product  (Figure  5).  The  calculated  spectra 
show  reasonable  agreement  with  the  observed  spectra,  indicating  this  model  may 
be  useful  to  explain  the  experimental  spectral  data. 


Predicted  Concentrations  of  Individual  Reaction 
Components 

The  concentration  curves  calculated  using  the  sequential  first-order  kinetic 
model  for  the  four  components  at  25  °C  for  the  first  3,000  s  of  the  TNT- 
hydroxide  reaction  are  shown  in  Figure  8.  TNT  concentration,  the  initial  com¬ 
ponent,  quickly  goes  to  extinction  as  indicated  by  the  thin  solid  line.  The  concen¬ 
tration  of  the  first  intermediate,  illustrated  by  the  dotted  line,  rises  and  falls 
quickly  in  the  reaction.  The  second  intermediate,  represented  by  the  bold  solid 
line,  grows  to  a  maximum  at  191  s  at  25  °C  and  very  slowly  decreases  over  time. 
The  final  concentration  of  the  product  slowly  increases  over  the  course  of  the 
reaction,  as  shown  by  the  dashed  line  on  Figure  8.  The  bottom  plot  illustrates  the 
concentration  curves  for  the  four  components  at  25  °C  for  the  first  500  s  of  the 
reaction  as  fit  by  the  sequential  first-order  model.  A  plot  of  the  entire  sampling 
time  was  not  illustrated,  as  the  smaller  features  were  obscured  in  the  plot. 


Other  Kinetic  Models 

In  this  work,  the  slow  formation  of  the  final  component  in  the  four- 
component  model  suggested  that  the  final  product  might  be  produced  in  a 
second-order  step.  A  reaction  scheme  for  this  model  is 

kl  A2  A? 

A  ->  B  ->  C  ->  D 

1st  order  1st  order  2nd  order 

for  which  -dCldt  =  k2.  [B]  -  k2  [C]2  and  the  other  terms  can  be  written  accord¬ 
ingly.  A  four-component,  first-order,  first-order,  second-order  kinetics  model  was 
applied  to  20  °C  spectral  data.  Figure  9  illustrates  the  spectra  of  the  four  com¬ 
ponents  of  the  TNT-hydroxide  reaction  using  this  model.  TNT,  first  intermediate, 
second  intermediate,  and  the  final  products  follow  very  similar  trends  at  20  °C 
for  this  model  at  25  °C  for  the  sequential  first-order  kinetics  model  (Figure  7). 
This  may  indicate  that  the  plotted  spectra  are  not  sensitive  to  the  models. 
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Figure  8.  Plot  of  concentration  data  for  the  four  components  at  25  °C  using 
sequential  first-order  kinetic  model 


The  kinetic  (concentration)  plots  for  20  °C  in  Figure  10  also  reflect  similar  trends 
as  those  at  25  °C  in  the  sequential  first-order  model  (Figure  8).  This  may  indicate 
that  both  models  are  working  reasonably  well  to  predict  the  spectral  and  kinetic 
data  of  the  TNT-hydroxide  reaction. 
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Figure  9.  Spectral  features  of  four  components  at  20  °C  using  first-order,  first- 
order,  second-order  kinetics  models 

Table  4  shows  the  rate  constants  that  were  calculated  using  both  four- 
component  models  that  were  described  previously.  The  values  for  k\  and  k3, 
calculated  using  the  sequential  first-order  kinetics  model,  have  equal  and  oppo¬ 
site  values.  The  third  rate  constant  is  negative  for  all  temperatures,  which  is  not 
physically  possible.  %2  were  in  the  order  of  10s,  and  SSE  was  calculated  using  the 
entire  experimental  matrix.  The  first-order,  first-order,  second-order  kinetics 
model  did  show  an  improvement  from  the  sequential  first-order  model  in  that  all 
the  rate  constants  were  positive.  Second-rate  constants  are  at  least  an  order  of 
magnitude  smaller  in  this  model.  %2  and  SSE  are  in  the  same  order  of  magnitude 
in  both  models.  This  table  shows  that  the  rate  constant  for  the  initial  step  in  the 
reaction  (ky)  is  much  larger  (faster)  than  the  second  rate  constant  (k2).  This 
corresponds  to  results  from  a  kinetics  study  that  followed  TNT  decay  using 
HPLC  (Felt,  Larson,  and  Hansen  2001a).  Rate  constants  were  calculated  in  that 
study  using  a  second  exponential  decay  to  model  TNT  transformation. 
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Table  4 

Rate  Constants  Calculated  Using  the  Sequential  First-Order,  Four- 
Component  Model  and  the  First-Order,  First-Order,  Second-Order, 
Four-Component  Model  (12'  represents  12-°C  data  that  had  been 
recalculated  using  the  15-°C  data  as  a  starting  point  for  the 
calculations) 

Sequential  First-Order  Model 

12 

15 

25 

ki 

0.0205 

0.0206 

0.0206 

0.0206 

k2 

0.0029 

0.00684 

0.04529 

0.04396 

ks 

-0.0204 

-0.0204 

-0.0205 

-0.0205 

X2 

9.82E+05 

6.60E+05 

8.84E+05 

6.78E+05 

SSE 

0.826 

0.556 

0.678 

0.571 

1-1-2  Model 

Temperature  (°C) 

12’ 

12 

15 

Hi 

ki 

0.02173 

0.1234 

0.0208 

0.0526 

0.076 

k2 

0.000116 

0.000171 

0.000143 

0.0002 

0.000234 

k3 

0.00676 

0.00855 

0.00583 

0.0228 

0.03151 

X2 

5.09E+05 

4.95E+05 

3.53E+05 

6.22E+05 

5.22E+05 

SSE 

0.429 

0.417 

0.297 

0.524 

0.465 

The  fits  for  the  four-component  models  are  interesting  but  not  perfected.  %2 
are  in  the  order  of  1 05,  were  reduced  from  1012,  and  the  values  were  set  at  the 
initial  step  of  the  calculation.  The  overall  fit  to  the  absorbance  differences  are  on 
the  order  of  0.04  AU,  and  most  below  0.008  AU,  but  these  values  are  above  the 
estimated  error  of  an  absorbance  reading  (0.002  AU).  Judging  from  the  X 
criteria,  the  exact  rate  model  is  as  yet  not  sufficiently  characterized.  The  exact 
rate  model  does  not  alter  the  basic  features  of  the  rise  and  fall  of  the  first  three 
components  and  the  growth  of  the  terminal  component  to  represent  the  spectra  at 
long  reaction  times.  The  four-component  models  based  on  PCA  therefore 
reasonably  map  out  the  course  of  the  reaction. 

A  kinetics  model  that  represents  multiple  products  may  be  illustrated  as: 

d  “►  e 

X 

a  ->  b  ->  c  "►  f  (23) 

* 

(TNT)  (Int.  1)  (Int.2)  g  h 

Other  kinetics  models  for  the  TNT-hydroxide  reaction  are  possible.  The 
observed  spectra  indicate  that  TNT  degrades  to  form  a  first  intermediate  (b)  that 
rapidly  degrades  to  form  a  second  intermediate  (c).  The  second  intermediate  may 
degrade  to  form  more  than  one  product,  which  was  indicated  by  a  previous  study 
that  was  conducted  in  our  laboratory.  Using  gel  permeation  chromatography 
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(GPC),  OH-TNT  final  reaction  products  were  separated  into  molecular  weight 
fractions  (Felt,  Larson,  and  Hansen  2001b).  TNT  final  products  varied  in 
molecular  weight  from  3,000  to  6,000  to  <100  Da.  GPC  results  do  not  contradict 
the  spectral  model  offered  here,  since  the  PCA  would  not  necessarily  differenti¬ 
ate  between  similar  chemical  compounds. 

It  should  be  noted  that  the  order  of  the  OH-TNT  reaction  has  not  been 
established  experimentally.  Knowing  the  order  of  the  reaction  could  verify  a 
kinetics  model  for  the  OH-TNT  reaction.  The  smaller  principal  components  of 
the  OH-TNT  reaction  could  be  identified  using  an  LC-MS  technique  after 
filtration  of  the  polymers.  The  overall  order  of  the  reaction  could  be  determined 
after  the  components  were  identified  and  the  order  for  each  component  was 
determined.  Kinetic  models  could  be  refined  and  tested,  if  the  order  of  the 
reaction  were  known.  Establishing  the  order  of  the  reaction  could  be  a  topic  for 
future  work. 
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Conclusions 


OH-TNT  reaction  mixtures  were  analyzed  over  time  using  UV/VIS  spectro¬ 
scopy  at  different  temperatures.  The  spectra  indicated  that  the  OH-TNT  reaction 
requires  48  h  to  come  to  completion  at  room  temperature,  long  after  TNT  itself 
has  degraded  (<  40  min).  An  intermediate  of  the  reaction  produced  a  dominant 
spectral  feature  with  a  maximum  of  0.8  AU  at  450  nm,  regardless  of  temperature. 
This  suggested  that  the  component  was  a  real  reaction  component,  and  the  spec¬ 
trum  containing  a  maximum  450-nm  feature  was  its  true  spectrum.  Results  indi¬ 
cate  an  initial  reaction  intermediate  is  rapidly  formed  after  the  addition  of  the 
base  that  quickly  transforms  into  a  polar,  secondary  intermediate.  The  secondary 
intermediate,  in  turn,  slowly  degrades  to  form  polar  final  products.  These  results 
were  supported  by  a  previous  kinetic  study  using  HPLC. 

The  observed  spectral  data  were  analyzed  by  FA.  PCA  of  the  observed  data 
indicated  that  there  were  six  principal  components  in  the  reaction,  four  of  which 
explained  most  of  the  variance  in  the  data  matrix.  A  sequential  first-order  model 
and  a  first-order,  first-order,  second-order  model  were  tried  and  were  judged  to 
be  interesting.  The  models  and  rate  measurements  were  abstract  and  exhibited 
trends  similar  to  rate  measurements  calculated  for  HPLC  data  by  a  second  expo- 
ential  decay  fit.  The  rate  order  of  the  OH-TNT  reaction  has  not  been  established 
experimentally.  Identification  of  individual  reaction  components  and  the  order  of 
the  reaction  would  help  determine  an  accurate  kinetics  model  for  the  OH-TNT 
reaction. 
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Appendix  A  Principal  Component  Analysis  (PCA)  for  Spectral  Data 


A1 


! - PRINCIPAL  COMPONENT  ANALYSIS  FOR  SPECTRA 

implicit  real*8(A-H,0-Z) 

CHARACTER*  1  TITL(80) 

DIMENSION  D(  1 00, 1 00), T(  1 00), WAV(  1 00) 

DIMEN  SION  DT(  1 00, 1 00),CO  V  AR(  1 00, 1 00) 

DIMENSION  ANORM(  1 00),DNORM ( 1 00, 1 00), Z(  1 00, 1 00) 

DIMENSION  DNORMT(  1 00, 1 00),EIGEN(  1 00) 

DIMENSION  V ( 1 00), VE(  1 00), VT ( 1 00) 

DIMENSION  C(100, 100) 

DIMENSION  ADIFF(  1 00, 1 00),DDIFF(  1 00, 1 00) 

DIMENSION  RE(  1 00),EMBE(  1 00) 

DIMENSION  CPRIME(  1 00, 1 00),DPRIME(  1 00, 1 00),ROW(  1 00, 1 00) 

DIMENSION  CPRIMET(  1 00, 1 00) 

! 

!  READ  DATA  MATRIX 

!  FIRST  READ  NUMBER  OF  ROWS  (IR)  THEN  NUMBER  OF  COLUMN  (IC) 

!  FOR  ABSORBANCE  VS.  TIME:  (OTHER  2D  VARIABLE  JUST  AS  ACCEPTABLE) 
!  THE  ROW  ENTRIES  ARE  THE  ABSORBANCES  AT  VARIOUS  TIMES 
!  AND  SUCCESSIVE  ROWS  ARE  ADDITIONAL  ABSORBANCES  AT  VARIOUS 
TIMES 

!  THE  TIME  DATA  SHOULD  BE  IN  THE  LAST  ROW  RECORD 
OPEN(3, FILE- HPDATA.DAT', STATUS-UNKNOWN') 

!  OPEN  OUTPUT  FILE,  TOO 

OPEN(10,FILE='PCA.OUT', STATUS- UNKNOWN') 

READ(3 , 1 0)IR,IC 
10  FORMAT(2I5) 

DO  1=1, IR 

READ(3,*)  (D(I,J),J=1,IC) 

ENDDO 

READ(3,*)  (T(J),J=1,IC) 

WRITE(6,25)IR,IC 
WRITE(  1 0,25)IR,IC 

25  FORMAT( IX, 'Principal  Component  Analysis 

/IX, 15,'  x',I5,'  DATA  MATRIX  READ  FROM  "HPDATA.DAT"') 

DO  1=1, IR 

WRITE(  1 0,30)(D(I,J),J=1,IC) 

30  FORM AT(  1  X,6(E  1 2.4, 1 X)) 

ENDDO 

WRITE(  10,31)  (T(J),J=1  ,IC) 

3 1  FORMAT(  1 X, 'TIMES  :’,/l  X,6(E  1 2.4, 1 X)) 

WRITE(10,*) 

!  READ  WAVELENGTH 
READ(3,*)(WAV(I),I=1  ,IR) 

READ(3,33)(TITL(I),I=1 ,80) 

33  FORMAT(80A1) 

CLOSE(3) 

!  COMPUTE  THE  TRANSPOSE  OF  THE  DATA  MATRIX; 

!  THEN  THE  COVARIANCE  MATRIX  DT.D  DIMENSION  NOW  ICxIR 
CALL  TRANSPOS(IR,IC,D,DT) 

CALL  M  ATMULT  (IC,IR,IC,DT,D,COV  AR) 
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WRJTE(10,29) 

29  FORM AT( IX, 'TRANSPOSED  DATA  MATRIX:') 

DO  1=1,  IC 

WRITE(  1 0,30)(DT  (I,J),  J=1  ,IR) 

ENDDO 
WRITE(  10,32) 

32  F0RMAT(1X, 'COVARIANCE  MATRIX:') 

DO  1=1  ,IC 

WRITE(1 0,30)(COVAR(I,J),J=1  ,IC) 

ENDDO 

!  NOW  NORMALIZE  THE  COVARIANCE  MATRIX 
!  FIRST  GET  THE  NORMALIZATION  CONSTANTS  FOR  EACH  COLUMN 
DO  1=1,  IR 
DO  J=1,IC 

ANORM(J)=ANORM(J)+(D(I,J)**2) 

ENDDO 

ENDDO 

!  WRITE(10,*)ANORM(1),ANORM(2),ANORM(3),ANORM(4) 

DO  1=1, IC 

ANORM(I)=  1  /ANORM(I) 

ANORM(I)=DSQRT(ANORM(I)) 

ENDDO 

!  WRITE(10,*)ANORM(1),ANORM(2),ANORM(3) 

!  COMPUTE  NORMALIZED  DATA  MATRIX  DNORM 
DO  1=1, IR 
DO  J=1,IC 

DNORM(I,J)=D(I,J)*ANORM(J) 

ENDDO 
ENDDO 
WRITE(  10,51) 

51  FORMAT('NORMALIZED  DATA  MATRIX:') 

DO  1=1, IR 

WRITE(10,30)(DNORM(I,J),J=1,IC) 

ENDDO 

!  AND  CALCULATE  THE  CORRELATION  MATRIX,  Z=DNORM*(DNORM)T 
CALL  TRANSPOS(IR,IC, DNORM, DNORMT) 

CALL  MATMULT(IC,IR,IC, DNORMT, DNORM, Z) 

WRITE(10,59) 

59  F0RMAT(1X, 'NORMALIZED  CORRELATION  MATRIX  AS  COR  TRIANGLE:') 
DO  1=1, IC 

WRITE(10,60)(Z(I,J),J=1 ,1) 

60  F0RMAT(13(F5.3,1X)) 

ENDDO 

i 

!  NOW  DETERMINE  THE  EIGENVECTORS  IN  THE  SPACE;&  THE  EIGENVALUES 
!  BY  ITERATION,  WE  SATISFY  [Z]*C1=EIG1*C1 
!  INTITIAL  ESTIMATES  OF  EIGENVECTOR 

WRITE(  1 0,*)"DETERMINE  EIGENVECTORS  BY  ITERATION" 

DO  1=1, IC 
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DO  J=1,IC 
C(I,J)=0.58 
ENDDO 
ENDDO 

!  BEGIN  ITERATION 
DEL=1.0D-9 
START=-1.0D-10 

!  COMPUTE  IC  EIGENVALUES  AND  EIGENVECTORS 
DO  121  N=1,IC 
DO  110  L=l,50 
DO  1=1, IC 
V(I)=C(N,I) 

ENDDO 

CALL  MATVEC(IC,IC,COVAR,V,VE) 

!  WRITE(  1 0, 1 04)(VE(J),J=1  ,IC) 

104  FORMAT(1X,10(1X,E10.5)) 

!  MORMALIZE  THE  VECTOR 

XNORM=O.ODO 
DO  1=1, IC 

XNORM=XNORM+(VE(I)**2) 

ENDDO 

XNORM=DSQRT(XNORM) 

DIFF=DABS(XNORM-START) 

WRITE(  1 0, 1 09)L,XNORM,DIFF 

109  FORM AT( IX, 'ITER#  ',13,'EST.  EIGENVALUE  =  ',F12.3,& 
'  CHANGE  =  ', FI 0.5) 

DO  1=1, IC 

VE(I)= VE(I)/XN  ORM 
ENDDO 

!  WRITE(  1 0, 1 04)(VE(J),  J=  1  ,IC) 

DO  1=1, IC 
C(N,I)=VE(I) 

ENDDO 

ST  ART=XN  ORM 
IF(DIFF.LT.DEL)GOTO  1 1 1 

110  CONTINUE 

111  CONTINUE 

WRITE(  1 0, 1 05)N,XNORM,(C(N,J),J=  1  ,IC) 

105  FORMAT( IX, 'EIGEN VALUE  ’,12,’  IS  ’,F12.5,/1X,& 
'VECTOR  IS  ’,3(1  X,F1 2.5)) 

!  END  ITERATION 

!  COMPUTE  RESIDUAL  MARIX  R,  =COVAR  -  LAM*C*C' 
EIGEN  (N)=XNORM 
DO  1=1, IC 

VT (I)= VE(I)  *  EIGEN (N) 

ENDDO 

CALL  VECMULT(IC,  1  ,VT, VE,ADIFF) 

!  DO  1=1, IC 

!  WRITE  (10,*)  (ADIFF(I,J),J=1  ,IC) 

!  ENDDO 
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!  SUBTRACT  FROM  COVARIANCE  MATRIX 
DO  1=1, IC 
DO  J=1,IC 

CO  V  AR(I,  J)=CO  V  AR(I ,  J)- ADIFF  (I,  J) 

ENDDO 

ENDDO 

121  CONTINUE 

!  PRINT  THE  COVARIANCE  MATIX  AFTER  ALL  IS  DONE  (SHOULD  BE  ZEROS) 
WRITE(10,*)"FINAL  COVARIANCE  MATRIX:" 

DO  1=1, IC 

WRITE(  10,1 22)(CO  V  AR(I,  J),  J=  1  ,IC) 

122  FORMAT(1X,10(F6.4,1X)) 

ENDDO 

i 

IVEC=IC 
IMPTIVEC=IC 
CUT=1.0D0/IC 
DO  1=1, IC 

IF  (EIGEN  (I)  .LT.  CUT)IMPTI  VEC=IMPTI  VEC- 1 
ENDDO 

!  OFFER  THE  OPPORTUNITY  TO  MODIFY  THE  NUMBER  OF  IMPORTANT 
EIGENVECTORS 

!  WRITE(6,*)'ENTER  NUMBER  OF  EIGENVECTORS  TO  OVERRIDE  CALC.  (II) 
0=KEEP 

!  READ(5 , 1 29)ICHOSE 
!  129  FORMAT(Il) 

!  IF(ICHOSE.EQ.O)GOTO  130 
!  IMPTIVEC=ICHOSE 

!  CALCULATE  REAL  AND  IMBEDDED  ERRORS 

!  RE  =  THAT  ACCOUNTED  FOR  BY  THE  OTHER  FACTORS  THAN  THOSE 
EXPLAINING  THE 
!  VARIANCE 
130  DO  I=1,IVEC 
K=I+1 

ANUM=0.0D0 
DO  J=K,IVEC 
ANUM=ANUM+EIGEN  ( J) 

ENDDO 

!  WRITE(6,*)I,J,ANUM 
ARG=ANUM/(IC*(IVEC-I)) 

RE(I)=DSQRT(ARG) 

ANUM=I*ANUM 

ARG=ANUM/(IR*IC*(IVEC-I)) 

EMBE(I)=DSQRT(ARG) 

ENDDO 

!  SOME  SCREEN  OUTPUT;  EIGENVALUES  AND  EIGENVECTORS 
WRITE(6, 1 3  5)IR,IC 

135  F0RMAT(1X,'DATA  MATRIX  HAS  ',13,'  ROWS  AND  ',13,'  COLUMNS') 
WRITE(6,1 39)IMPTIVEC 

139  F ORM AT ( 1  X,'PC A  EXTRACTED  ',13,'  IMPORTANT  EIGENVECTORS') 
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DO  I=1,IVEC 
WRITE(6, 143)1,  EIGEN(I) 

WRITE(6, 1 44)(C(I,J),J=1,IC) 

WRITE(6, 1 45)RE(I),EMBE(I) 

WRITE(  10,1 43)1, EIGEN(I) 

WRITE(  1 0, 1 44)(C(I,J),J=  1  ,IC) 

WRITE(  1 0, 1 45)RE(I),EMBE(I) 

143  FORM AT(1X,'EIGEN VALUE  #  ',I2,5X,F12.3) 

144  F0RMAT(1X,  10(F7.2,1X)) 

145  F0RMAT(1X,'REAL  ERROR=  ',F6.2,'  EMBEDDED  ERROR=  ',F6.2) 

ENDDO 

! 

!  COMPUTE  THE  DATA  MATIX  FROM  THE  PCA  SOLUTION 
!  D'=C*R  USING  THE  IMPTIVEC  EIGENVECTORS  AND  EIGENVALUES 
!  CONSTRUCT  THE  EIGENVECTOR  APPROXIMATION 
DO  1=1, IMPTIVEC 
DO  J=1,IC 
CPRIME(I,J)=C(I,J) 

ENDDO 

ENDDO 

!  WRITE(10,*)"CPRIME" 

!  DO  1=1, IMPTIVEC 

!  WRITE(  10,1 44)(CPRIME(I,J),J=  1  ,IC) 

!  ENDDO 

!  AND  GENEERATE  THE  ROW  MATRIX  FROM  THE  IMPORTANT 
EIGENVECTORS 
!  ACCORDING  TO  R=d*Ct 

CALL  TRANSP0S(IMPT1VEC,IC,CPRIME,CPRIMET) 

!  WRITE(  1 0,  *)"CPRIME  TRANSPOSE" 

!  DO  1=1, IC 

!  WRITE(10,144)(CPRIMET(I,J),J=1, IMPTIVEC) 

!  ENDDO 

!  AND  MULTIPLY  BY  THE  EIGENVECTORS  TO  GENERATE  THE  ROW  MATRIX 
CALL  MATMULT(IR,IC, IMPTIVEC, D,CPRIMET, ROW) 

WRITE(10, ^"RECONSTRUCTED  ROW  MATRIX  BASED  ON  IMPORTANT 
FACTORS" 

WRITE(10,146) 

146  FORM AT(  IX, 'LAMBDA  1  2  3  4  5  6') 

DO  1=1, IR 

WRITE(  10,1 47)  WA  V(I),(RO  W(I,  J),  J=  1 , IMPTIVEC) 

ENDDO 

1 47  FORMAT(1X,F4.0,3X,  1 0(F7.2, 1 X)) 

!  AND  FINALLY,  THE  RECONSTRUCTED  DATA  IS  ROW*CPRIME 
CALL  MATMULT(IR, IMPTIVEC, IC, ROW, CPRIME,DPRIME) 

WRITE(  1 0,*)"RECONSTR.  DATA  MATRIX  BASED  ON  ", IMPTIVEC," 
IMPORTANT  FACTORS" 

DO  1=1, IR 

WRITE(  1 0,30)(DPRIME(I,J),J=  1  ,IC) 

ENDDO 

!  AND  THE  DIFFERENCE  BETWEEN  THE  DATA  AND  ITS  RECONSTRUCTION 
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DO  1=1, IR 
DO  J=1,IC 

DDIFF  (I,  J)=D(I,  J)-DPRIME(I,  J) 

ENDDO 

ENDDO 

WRITE(10,*)"DIFF.  BETW.  DATA  AND  RECONSTRUCTION  =  ERROR" 

DO  1=1, IR 

WRITE(1 0, 1 5 1  )(DDIFF(I,J),J=1  ,IC) 

151  FORMAT(1X,10(F6.3,1X)) 

ENDDO 

i 

!  ISSUE  OUTPUT  FOR  READING  TO  TEST  DIRECTED  VECTORS 
OPEN(l  1, FILE- TESTVEC.DAT', STATUS-UNKNOWN') 

WRITE(  11,1 0)IR,IC 
!  10  FORMAT(2I5) 

DO  1=1, IR 

WRITE(1 1,*)(D(I,J),J=1,IC) 

ENDDO 
DO  1=1, IR 

WRITE(1 1,*)(DPRIME(I,J),J=1,IC) 

ENDDO 

WRITE(1 1,10)IVEC,IMPTIVEC 
DO  I=1,IVEC 
WRITE(1 1  ,*)EIGEN(I) 

ENDDO 
DO  I=1,IVEC 

WRITE(1 1,*)(C(I,J),J=1,IC) 

ENDDO 
DO  1=1, IR 

WRITE(1 1  ,*)(ROW(I,J),J=l  ,IMPTIVEC) 

ENDDO 

WRITE(1 1,*)(WAV(I),I=1,IR) 

WRITE(1 1,*)(T(I),I=1,IC) 

WRITE(10,*)"DATA  FOR  MODEL  TESTING  WRITTEN  TO  'TESTVEC.DAT'" 
WRITE(10,233)(TITL(I),I=1,80) 

WRITE(1 1,233)(TITL(I),I=1,80) 

233  FORMAT(1X,80A1) 

CLOSE(ll) 

WRITE(6,*)"OUTPUT  WRITTEN  TO  ’PCA.OUT'" 

STOP 

END 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

! 

!  SUBROUTINES 
! 

!  COMPUTE  THE  TRANSPOSE  OF  A  RECTANGULAR  MATRIX  !!!!!!!!!!!!!!!!!!!!!!!!!! 
SUBROUTINE  TRANSPOS(II,JJ,A,B) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00, 1 00),B(  1 00, 1 00) 

!  WRITE(6,*)"TRANSPOSE",II,JJ 
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DO  1=1,11 
DO  J=1,JJ 
B(J,I)=A(I,J) 

!  WRITE(10,*)B(J,I) 

ENDDO 
ENDDO 
RETURN 
END 


!  COMPUTE  THE  MATRIX  PRODUCT  OF  AN  IxK  AND  KxJ  ARRAYS  TO  PRODUCE 
A 

!  THE  GENERAL  PRODUCT  MATRIX  OF  DIMENSION  IxJ 
SUBROUTINE  MATMULT(II,KK,JJ, A, B,AB) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(1 00, 1 00), B(  1 00, 1 00),AB(1 00, 1 00) 

!  BY  THE  ROWxCOLUMN  METHOD;  AN  II  BY  KK  MATRIX  IS  MULTIPLIED  BY  A 
!  KK  BY  JJ  PRODUCING  AN  II  BY  JJ  MATRIX 
DO  1=1,11 
DO  J=1,JJ 
SUM=0.0D0 
DO  K=1,KK 

SUM=SUM+(A(I,K)*B(K,J)) 

!  WRITE(  1 0,40)1,  J,K,B(I,K),A(K,J), SUM 

!  40  F0RMAT(1X,3I5,1X,3(F12.4,1X)) 

ENDDO 

AB(I,J)=SUM 

ENDDO 

ENDDO 

RETURN 

END 

! 

!  MULTIPLY  A  VECTOR  (B)  BY  A  MATRIX  (A)  [DOT  PRODUCT] 

!  COLUMNS  OF  VECTOR  MUST  EQUAL  ROWS  OF  MATRIX  (IR) 

!  OUTPUT  VECTOR  (C)  HAS  IR  ROWS 
SUBROUTINE  MATVEC(IR,IC,A,B,C) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00, 1 00), B(  1 00),C(1 00) 

DO  1=1,  IC 
SUM=0.0D0 
DO  J=1,IR 

SUM=SUM+(A(I,J)*B(J)) 

!  WRITE(  1 0, 30)1,  J,A(I,J),B(J), SUM 
!  30  F0RMAT(1X,I5,I5,3(F12.5,1X)) 

ENDDO 

C(I)=SUM 

ENDDO 

RETURN 

END 
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!  MULTIPLY  TWO  VECTORS  (A)  AND  (B);  (RxC)*(CxR)  =  (RxR) 
!  TO  PRODUCE  A  MATRIX 

SUBROUTINE  VECMULT(II,JJ,A,B,C) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00), B(  1 00), C(  1 00, 1 00) 

DO  1=1,11 
DO  J=1,II 
C(I,J)=A(I)*B(J) 

!  WRITE(10,30)I,J,A(I),B(J),C(I,J) 

!  30  F0RMAT(1X,2I5,1X,3F12.5) 

ENDDO 

ENDDO 

RETURN 

END 
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!  PROGRAM  TO  TEST  VECTORS  FROM  PCA 

!  VERSION  FOR  TWO  FIRST  ORDER  FOLLOWED  BY  SECOND  ORDER  KINETICS 
!  TREATMENT 

!  MODIFIED  LEAST-SQUARES  PROCEDURE 
!  LOCAL  NAME  IS  TESTVEC  VERSION  7 
! 

IMPLICIT  REAL*8(A-H,0-Z) 

CHARACTER*  1  TITLE 

DIMEN  SION  D(  1 00, 1 00),DPRIME(  1 00, 1 00),  C(  1 00, 1 00),CBAR(  1 00, 1 00) 
DIMENSION  DD(1 00, 1 00),ROWT(  1 00, 1 00),WAV(  1 00) 

DIMEN  SION  EIGEN(  1 00),ALAM(  1 00, 1 00),ALAM  IN  V(  1 00, 1 00),ROW(  1 00, 1 00) 
DIMENSION  TT(  1 00, 1 00),TL(  1 00),XFORM(  1 00, 1 00),TIN V(  1 00, 1 00) 

DIMENSION  DATA(  1 00),UNITY(  1 00),UNIQUE(  1 00),SELF(  1 00),TEST(  1 00) 
DIMENSION  TITLE(80) 

DIMENSION  A(  1 00,1 00), A A(  1 00) 

DIMENSION  SCRATCH(500) 

DIMENSION  TIME(  1 00),CONC(  1 0, 1 00) 

i 

DIMENSION  G(  1 00, 1 0) 

DIMENSION  V(  1 00),RATEK(  1 0),CN(  1 0),TAU(  1 0, 1 0),ZETA(  1 0, 1 0) 

DIMENSION  ACALC(  1 00, 1 00),Q(  1 00, 1 00) 

DIMENSION  SIG(  1 00, 1 00),D2XIDADB(  1 00, 1 00),DXIDV(  1 00) 

DIMENSION  ADEL(  100,1 00),DELA(  1 00) 

DIMENSION  DADVAR(  100, 100,50) 

DIMENSION  BLO(  1 0),BHI(  1 0),SCALE(  1 0) 

!  RETRIEVE  DATA  FROM  UNIT  1 1  TESTVEC.DAT',  MADE  BY  PCA.F90 
OPEN(  1 1 , FILE- TESTVEC. DAT',STATUS='UNKNOWN') 

!  OUTPUT  IS  UNIT  10 

OPEN(  1 0,FILE='TESTVEC.OUT', STATUS- UNKNOWN') 

READ(1 1,10)IR,IC 
10  FORMAT(2I5) 

WRITE(6,*)IR,IC 

WRITE(10,*)IR,'  ROWS', IC,'  COLUMNS' 

DO  1=1, IR 

READ(  1 1  ,*)(D(I,J),J=1  ,IC) 

ENDDO 
!  GET  LIMITS 
ALOW=1.0D+25 
AHI=1.0D-25 
DO  1=1, IR 
DO  J=1,IC 

IF(D(I,J).GT.AHI)AHI=D(I,J) 

IF(D(I,J).LT.ALOW)ALOW=D(I,J) 

ENDDO 

ENDDO 

WRITE(  1 0, 1 4)ALOW,AHI 

14  FORMAT( IX, 'ABSORBANCE  RANGE  FROM  ',F7.3,'  TO  ’,F7.3) 

!  WRITE(  1 0,*)'DATA  MATRIX' 

!  DO  1=1, IR 
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!  WRITE(10,1 1)(D(I,J),J=1,IC) 

!  ENDDO 
DO  1=1,  IR 

READ(1 1,*)(DPRIME(1,J),J=1,IC) 

ENDDO 

!  WRITE(10, ^'RECONSTRUCTED  DATA  MATRIX' 

!  DO  1=1, IR 

!  WRITE(10,1 1)(DPRIME(I,J),J=1,IC) 

!  ENDDO 

11  FORMAT(1X,6(F10.6,1X)) 

READ(  11,1 0)IVEC,IMPTI  VEC 
WRITE(  1 0,8)IVEC,IMPTIVEC 
WRITE(6,8)IVEC,IMPTIVEC 

8  F0RMAT(1X,I3,'  VECTORS  WITH  ',13,'  IMPORTANT  ONES') 

DO  I=1,IVEC 
READ(1 1,*)EIGEN(I) 

ENDDO 

WRITE(1 0,*)'EIGENVALUES' 

DO  1=1, 1  VEC 
WRITE(  1 0,*)EIGEN(I) 

!  WRITE(6,  *)EIGEN  (I) 

ENDDO 
DO  I=1,IVEC 
READ(1 1  ,*)(C(I,J),J=1,IC) 

ENDDO 

WRITE(10,*)'EIGENVECTORS  FOR  IMPORTANT  ONES' 

DO  I=1,IMPTIVEC 
WRITE(  1 0,*)(C(I,J),J=1  ,IC) 

ENDDO 

!  READ  IMPTIVEC  RECONSTRUCTED  ROW  VECTORS  (FOR  TESTS) 

DO  1=1, IR 

READ(1 1  ,*)(ROW(I,J),J=l , IMPTIVEC) 

ENDDO 

!  READ  IN  THE  WAVELENGTHS  FOR  EACH  ROW  OF  THE  ROW  MATRIX 
READ(1 1  ,*)(WAV(I),I=1  ,IR) 

!  WRITE(6,*)(WAV(I),I=1,IR) 

!  READ  IN  THE  TIMES  FOR  THE  EIGENVECTORS 
READ(11,*)(TIME(I),I=1,IC) 

!  WRITE(6,*)(TIME(I),I=1  ,IC) 

READ(1 1,3)(TITLE(I),I=1,80) 

WRITE(6,3)(TITLE(I),I=1 ,80) 

3  FORMAT(80A1) 

CLOSE(ll) 

!  GENERATE  DATA'  MATRIX  JUST  TO  BE  SURE  ALL  IS  WELL 
CALL  M  ATMULT  (IR,IC  ,IR,RO  W,C  ,DD) 

DO  1=1, IR 
DO  J=1,IR 

DIFF=DD(I,J)-DPRIME(I,J) 

IF(DIFF.GT.  1 D-5 .OR.DIFF.LT.-l  D-5)  GO  TO  60 
ENDDO 
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ENDDO 
GO  TO  65 

60  WRITE(6,*)"ERROR  ON  DATA  READ" 

WRITE(6,*)DD(I,J),DPRIME(I,J),DIFF 
GO  TO  1001 

!  FILL  OUT  EIGENVALUE  DIAGONAL  MATRIX  (IVEC  x  IVEC) 

65  DO  1=1, IVEC 
DO  J=1,IVEC 
ALAM(I,J)=0.0D0 
ALAMINV(I,J)=0.0D0 
IF(I.EQ.J)ALAM(I,J)=EIGEN(I) 

IF(LEQ.J)ALAMINV(I,J)=1.0D0/EIGEN(I) 

ENDDO 

ENDDO 

WRITE(10,*)'IMPORTANT  EIGENVECTOR  LAMBDA  ARRAY' 

DO  I=1,IMPTIVEC 

WRITE(  1 0, 1 3)(ALAM(I,J),J=1  ,IMPTI  VEC) 

!  WRITE(6, 1 3)(ALAMINV(I,J),J=1  ,IVEC) 

13  FORMAT(1X,10(F8.3,1X)) 

ENDDO 

!  GENERATE  UNITY,  UNIQUENESS,  AND  A  SELF-ROW  VECTOR 
DO  1=1, IR 
UNITY(I)=1.0D0 
UNIQUE(I)=0.0D0 
SELF(I)=ROW(I,l) 

ENDDO 

!  TESTVECTOR  TT  =  (ALAM)-l  *  ROW  *  TESTVECTOR 
CALL  TRANSPOS(IR,IMPTIVEC,ROW,ROWT) 

CALL  MATMULT(IMPTIVEC,IMPTIVEC,IC,ALAMINV,ROWT,TT) 

! 

!  FIRST,  TEST  FIRST  ROW  VECTOR  FROM  THE  RECONSTRUCTED  [R]  OF  [R][C] 
WRITE(6,*)"DATA  (ROW  1)  TEST  =  SELF  TEST" 

WRITE(  1 0,*)"DATA  (ROW  1)  TEST  =  SELF  TEST" 

J=1 

DO  1=1, IR 
DATA(I)=ROW(I,J) 

ENDDO 

!  TL  IS  THE  TRANSFORMATION  FOR  THE  FIRST  COLUMN  OF  THE  ROW  MATRIX 
CALL  MATMULT(IMPTIVEC,IR,1,TT,DATA,TL) 

!  TEST  =  [ROW]*TL 

CALL  MATMULT(IR,IMPTI  VEC,  1 , ROW, TL, TEST) 

WRITE(6,23)  (ROW(I,J),I=l,IR) 

WRITE(6,*)'FIRST  COLUMN  OF  THE  RECONSTRUCTED  ROW  MATRIX:' 
WRITE(6,23)  (TEST(I),I=1,IR) 

WRITE(  10,23)  (ROW(I,J),I=l,IR) 

WRITE(  1 0,*)'FIRST  COLUMN  OF  THE  RECONSTRUCTED  ROW  MATRIX:' 
WRITE(  10,23)  (TEST(I),I=1,IR) 

23  FORMAT(l  1X,10(F5.2,1X)) 

300  CONTINUE 
!  UNITY 
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CALL  MATMULT  (IMPTI  VEC,IR,  1  ,TT, UNITY, TL) 

CALL  MATMULT(IR,IMPTIVEC,1 , ROW, TL, TEST) 

DIFF=0.0D0 
DO  J=1,IR 

DIFF=DIFF+(UNITY  (J)-TEST  (J))*  *2 
ENDDO 

WRITE(6,2 1  )DIFF 
WRITE(  1 0,2 1  )DIFF 

21  F0RMAT(1X, "UNITY  VECTOR:  (ALL  ONES)  SSE  =  ",F6.4) 

!  MAKE  UNITY,  SELF  AND  UNIQUENESS  CHECKS 

!  SELF 

DO  L=l, IMPTI VEC 
DO  1=1, IR 
SELF  (I)=RO  W  (I,L) 

ENDDO 

CALL  MATMULT  (IMPTI  VEC  ,IR,  1  ,TT,SELF,TL) 

CALL  MATMULT  (IR, IMPTI  VEC,  1 , ROW, TL, TEST) 

DIFF=0.0D0 
DO  J=1,IR 

DIFF=DIFF+(SELF(J)-TEST(J))**2 

ENDDO 

WRITE(6,29)L,DIFF 
WRITE(  1 0,29)L,DIFF 

29  F0RMAT(1X,"SELF  VECTOR:  ROW", 12,"  SSE  =  ”,F6.4) 

ENDDO 

!  UNIQUENESS  (CYCLE  OVER  ALL  EIGENVECTORS) 

DO  1=1,  IC 
UNIQUE(I)=1 .0D0 
K=I-1 

IF(I.GT.  1  )UNIQUE(K)=0.0D0 

CALL  MATMULT (IMPTIVEC,IR,  1  ,TT, UNIQUE, TL) 

CALL  MATMULT  (IR, IMPTI  VEC,  1 , ROW, TL, TEST) 

DIFF=0.0D0 
DO  J=1,IC 

DIFF=DIFF+(UNIQUE(J)-TEST  (J))*  *2 
ENDDO 

WRITE(6,22)I,DIFF 
WRITE(  1 0,22)I,DIFF 

22  F0RMAT(1X,I3,"TH  UNIQUE  VECTOR  SSE  =",F6.4) 

!  WRITE(6,23)  (TEST(J), J=  1  ,IR) 

!  WRITE(  10,22)  I,DIFF 
!  WRITE(10,23)  (TEST(J),J=1,IR) 

ENDDO 

!  RETRIEVE  SPECTROSCOPY  INFORMTION  »  INFO.DAT  CONTAINS 
!  A  TITLE 

!  TEST  VECTORS  1  ...  IMPTIVEC  AT  IR  WAVELENGTHS 
!  INSTRUCTIONS  (311)  ILSQ,  NCOM,  ICONTINUE 

!  THE  BALANCE  OF  THE  INFORMATION  HAS  COME  FROM  TESTVEC.DAT 
OPEN(8, FILE- INFO.DAT', STATUS-UNKNOWN') 
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READ(8,20 1  )(TITLE(I),I=  1 ,80) 

201  FORMAT(80A1) 

WRITE(6,*) 

WRITE(10,*) 

WRITE(6,20 1  )(TITLE(I),I= 1,80) 

WRITE(  1 0,20 1  )(TITLE(I),I=  1,80) 

!  READ  IN  IMPTIVEC  SPECTRA  (GUESSES),  AT  IR  WAVELENGTHS 
WRITE(6,223)(WAV(I),I=1  ,IR) 

223  FORMAT(l  1X,10(F4.0,2X)) 

DO  J=l, IMPTIVEC 
READ(8,*)(A(J,I),I=1  ,IR) 

ENDDO 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!  IF  ILSQ  =  1  PERFORM  FURTHER  ANALYSIS  FOR  THE  MODEL  CODED  HERE 
!  IF  NCOM  =  NUMBER  OF  COMPONENTS  (DEFAULT  =  IMPTIVEC) 

!  IF  ICONTINUE  >  0  DO  2A(ICONTINUE)  CYCLES 

!  IF  IREAD  =1,  READ  OLD  PARAMETERS  AND  BEGIN  THERE;  WRITE  NEW  AT  END 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

READ(8,288)ILSQ, NCOM, ICONTINUE, IREAD 
288  FORM AT(1 1,1 1,1 1,11) 

IF(NCOM.EQ.O)NCOM=IMPTIVEC 

ICY=2**(ICONTINUE) 

!!  TEST  THESE  ROW  MATRICES  -  BIG  LOOP 
DO  279  LL=1, IMPTIVEC 
DO  1=1, IR 
AA(I)=A(LL,I) 

ENDDO 

CALL  MATMULT(IMPTIVEC,IR,  1  ,TT,AA,TL) 

CALL  MATMULT(IR,IMPTIVEC,1 , ROW, TL, TEST) 

!  WRITE  COMPARISON 
DIFF=0.0D0 
DO  J=1,IR 

DIFF=DIFF+(AA(J)-TEST(J))**2 

ENDDO 

WRITE(6,33)LL,DIFF 
WRITE(6,26)  (ROW(I,LL),I=l  ,IR) 

WRITE(6,24)  (AA(I),I=1,IR) 

WRITE(6,25)  (TEST(I),I=1,IR) 

33  FORMAT(lX,I3,'TH  SPECTRAL  VECTOR  SSE=  ’,F7.5) 

WRITE(10,33)LL,DIFF 
WRITE(  10,26)  (ROW(I,LL),I=l,IR) 

WRITE(  10,24)  (AA(I),I=1,IR) 

WRITE(  10,25)  (TEST(J),J=1,IR) 

26  FORMAT(lX,'ROW  ’,15(F5.3,1X)) 

24  FORMAT(lX, ’GUESS ',15(F5.3, IX)) 

25  FORM AT( IX, TEST  ’,15(F5.3,1X)) 

!  NEW  IMPLIED  ROW  MATRIX  IS  AA 

!  SAVE  THE  TRANSFORM  (XFORM),  EACH  HAS  DIMENSION  IMPTIVEC  x  1 
DO  KK=1, IMPTIVEC 
XFORM(LL,KK)=TL(KK) 


B6 


Appendix  B  Program  to  Test  Vectors  from  PCA  Version 


ENDDO 

WRITE(  1 0,278)(XFORM(LL,KK),KK=1 , IMPTIVEC) 

!  END  LOOP  OVER  LL IMPTIVEC  VECTORS 
279  CONTINUE 
278  F0RMAT(1X,6(F12.5,1X)) 

DO  LL=1, IMPTIVEC 
DO  KK=1, IMPTIVEC 
TINV(LL,KK)=XFORM(LL,KK) 

ENDDO 

ENDDO 

CALL  MINV(TINV,IMPTIVEC, 200, SCRATCH, DET,  1  .OD-4, 1,1) 

!  WRITE(6,*)'IN  VERSE  TRANSFORM  MATRIX  (IMPTIVEC  x  IMPTIVEC)’ 

WRITE(  1 0,*)’INVERSE  TRANSFORM  MATRIX  (IMPTIVEC  x  IMPTIVEC)’ 

DO  LL=1, IMPTIVEC 

!  WRITE(6,27 8)(TIN  V  (LL,KK),KK=1 , IMPTIVEC) 

WRITE(  1 0,27 8)(TINV  (LL,KK),KK=  1 , IMPTIVEC) 

ENDDO 

!  AND  GENERATE  ROTATED  EIGENVECTORS 

CALL  MATMULT(IMPTIVEC, IMPTIVEC, IC,TINV,C,CBAR) 

WRITE(6,*)’OLD  EIGENVECTORS' 

WRITE(10,*)’OLD  EIGENVECTORS' 

DO  1=1, IMPTIVEC 
WRITE(6,282)I,(C(I,J),J=1  ,IC) 

WRITE(10, 282)1, (C(I,J),J=1,IC) 

282  F0RMAT(1X,I1,':,,15(F5.2)) 

ENDDO 

WRITE(6,*)'ROTATED  EIGENVECTORS' 

WRITE(10,*)'ROTATED  EIGENVECTORS' 

DO  1=1, IMPTIVEC 

WRITE(6, 282)1, (CBAR(I,J),J=1,IC) 

WRITE(  1 0, 282)1, (CBAR(I,J),J=  1  ,IC) 

ENDDO 

WRITE(6,*)'TIMES: ' 

WRITE(  1 0,  *)'TIMES : ' 

WRITE(6,283)(TIME(I),I=1,IC) 

WRITE(  1 0,283)(TIME(I),I=  1  ,IC) 

283  FORMAT(l X,  1 0(F7.0)) 

IF(ILSQ.EQ.0)GO  TO  1001 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!  A  TO  B  TO  C  ...  FITTING  ROUTINE;  UP  TO  NINE  IN  SEQUENCE 

!  FIND  THE  RATE  CONSTANTS  AND  THE  TRANSFORMATION  MATRIX  THAT  FITS  THE 
DATA 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

WRITE(6,*)'PERFORMING  NON-LINEAR  LEAST-SQUARES  ON  SEQENTIAL  MODEL' 

!  ESTIMATE  TAU  ELEMENTS  OF  TRANSFORMATION  MATRIX  FROM  TEST  VECTORS 
!  NUMBER  OF  VARIABLES  =  IMPTIVECA2  +  (IMPTIVEC- 1) 

K=1 

DO  1=1, IMPTIVEC 
DO  J=l, IMPTIVEC 
V(K)=TINV(I,J) 
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K=K+1 
ENDDO 
ENDDO 
NK=NCOM-l 
DO  1=1, NK 

!  DO  NOT  LET  INITIAL  GUESSES  OF  RATEK'S  BE  EQUAL 
V(K)=1.0D-2+(K*3.0D-3) 

K=K+1 

ENDDO 

!  WILD  GUESS  AT  RATE  CONSTANTS 
V(37)=0.0221 
V(38)=0.03 
V(39)=0.0455 
V(40)=0.1420 
V(41)=0.00000000001 

! 

NVAR=K-1 

WRITE(6,287)NVAR,NCOM,ICY 

287  F0RMAT(1X, 13,' VARIABLES  ',13,' COMPONENTS  ’,13,'  CYCLES') 

!  READ  OLD  PARAMETERS  IF  IREAD  =1 
IF(IREAD.EQ.  1  )READ(8,992)F 
IF(IREAD.EQ.1)READ(8,993)(V(I),I=1,NVAR) 

!  SET  DELV  FOR  DIFFERENTIALS  ON  ANY  VARIABLE 
DELV=1.0D-10 

!  PLACE  VARIABLES  INTO  TAU  AND  RATEK  ARRAYS 
DO  1=1  ,N  VAR 

CALL  VARFIG(I,NCOM, TAU, RATEK, V) 

ENDDO 

!  WRITE(  1 0,278)(RATEK(I),I= 1  ,NK) 

!  APPROXIMATE  SIGMAS  FOR  ABSORBANCES 
DO  1=1, IR 
DO  J=1,IC 
SIG(I,J)=4.0D-3 
ENDDO 
ENDDO 
NCY=0 

!  ENTER  ITERATION  ROUTINE . 

FUDGE=0. 1 D+003 
IF(IREAD.EQ.  1  )FUDGE=F 
CHISQ=1.0D+150 
IFLAG=0 

285  CHIOLD=CHISQ 
CHISQ=0.0D0 
NCY=NCY+1 

!  LOOP  OVER  ALL  VARIABLES 
DO  NV=1,NVAR 

!  CALC  Q  AND  XIA2;  LOOP  OVER  IR  ROWS(LAMBDAS)  AND  IC  COLUMNS(TIMES) 
DO  1=1, IR 
DO  J=1,IC 

!  THE  SECULAR  EQUATIONS  ARE  CODED  HERE 
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!  CALUCLATE  CONCENTRATIONS  OF  THE  SPECIES  (CN(1) ..  CN(NCOM)) 

CALL  CONCEN(NCOM,J,RATEK,TIME,CN,CONC) 

!  IF(J.EQ.6)WRITE(6,278)(CONC(K,J),K=l,NCOM) 

!  NOW  COMPUTE  THE  ZETA 

CALL  MODEL(I,NCOM,CN,TAU,ZETA, ROW, TERM) 

ACALC(I,J)=TERM 

Q(I,J)=D(I,J)-ACALC(I,J) 

!  NORMALIZE  THE  DATA 

CHISQ=CHISQ+(Q(I,J)/SIG(I,J))**2 
!  WRITE(6,*)CHISQ 
ENDDO 
ENDDO 
ENDDO 

WRITE(6,293)NCY,CHISQ,CHIOLD, FUDGE 
293  F0RMAT(1X,'CY#',I3,'CHISQ  =',E14.7,'  CHIOLD  =',E14.7,'  FUDGE  =',E14.7) 

WRITE(  1 0,293)NCY,CHISQ, CHIOLD, FUDGE 
IF(NCY.EQ.  1  )GOTO  299  !  FIRST  CYCLE 

!  DECISION  POINT:  IF  CHISQCHIOLD,  DECREASE  FUDGE,  ACCEPT  SHIFTS  AND  ITERATE 
!  IF  CHISQ>CHIOLD,  INCREASE  FUDGE,  REJECT  SHIFTS  AND  CYCLE 

!  IFLAG=1  FROM  PREV  CYCLE  WHERE  XIA2  INCREASED 
!  STOP  AFTER  NCY=30  IF  CHISQ<CHIOLD  OR  IF  CHISQ  <  ID- 10 
IF(ICY.EQ.l)GOTO  1001 
IF(IFLAG.EQ.O.AND.NCY.GE.ICY)GOTO  415 
IF(CHISQ.LT.  1  .OD-0 1 0)GOTO  4 1 5 
IF (IFLAG  .EQ .  1  )GOTO  299 
IF(CHISQ.GT.CHIOLD)GOTO  297 
FUDGE=FUDGE/1 .25D0 
GOTO  299 

297  FUDGE=FUDGE*  1 5 .0D0 
IFLAG=1 
DO  NV=1,NVAR 

V  (NV)=V(N  V)-DELA(NV) 

CALL  VARFIG(NV,NCOM,TAU,RATEK,V) 

ENDDO 
GOTO  285 

!  INCREMENT  A  VARIABLE  AND  SET  THE  TAU  AND  RATEK  VALUES 
299  IFLAG=0 
!  ENTER  LEAST-SQUARES 
DO  NV=1  ,NVAR 

V  (N  V)=V  (N  V)+DELV 

CALL  VARFIG(NV,NCOM,TAU,RATEK,V) 

!  CALCULATE  THE  DERIVATIVE  DY/DVAR 
DXIDV  (NV)=0.0D0 
DO  1=1  ,IR 
DO  J=1,IC 

CALL  CONCEN(NCOM,J, RATEK, TIME, CN,CONC) 

CALL  MODEL(I,NCOM,CN, TAU, ZETA, ROW, TERM) 

ADEL(I,J)=TERM 

DADVAR(I,J,NV)=((ACALC(I,J)-ADEL(I,J))/DELV) 

DXIDV(NV)=DXIDV(NV)+Q(I,J)*(1.0D0/(SIG(I,J))**2)*DADVAR(I,J,NV) 
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ENDDO 

ENDDO 

!  CALC  DIFF  OF  XI  WRT  VARIABLE  (A  VECTOR  DXIDV) 

!  WRITE(6,*)DXIDV(NV) 

DXIDV(NV)=-2.0D0*DXIDV(NV) 

V  (N  V)=V  (N  V)-DELV 

CALL  VARFIG(NV,NCOM,TAU,RATEK,V) 

ENDDO 

!  CALC  APPROX  SECOND  DERIV  DD  XI  WRT  VAR  A  &  VAR  B  (A  MATRIX  D2XIDADB) 
!  THIS  IS  APPROXIMATELY  THE  HESSIAN  OR  CURVATURE  MATRIX 
DO  NV=1,NVAR 
DO  NVB=1  ,NVAR 
D2XIDADB(NV,NVB)=0.0D0 
DO  1=1, IR 
DO  J=1,IC 

ARG=(1.0D0/(SIG(I,J))**2)*DADVAR(I,J,NV)*DADVAR(I,J,NVB) 

D2XIDADB(NV,NVB)=D2XIDADB(NV,NVB)+ARG 

ENDDO 

ENDDO 

D2XIDADB(NV,NVB)=2.0D0*D2XIDADB(NV,NVB) 

ENDDO 

ENDDO 

!  AND  SET  UP  LEVINSON-MARQUARDT  FUDGE  METHOD 
DO  I=1,NVAR 
DO  J=1,NVAR 
IF(I.NE.J)GOTO  295 

D2XIDADB(I,J)=D2XIDADB(I,J)*(  1  .ODO+FUDGE) 

295  ENDDO 
ENDDO 

!  SET  UP  MATRICES  FOR  GUASSIAN-JORDAN  ELIMINATION  SOLUTION  OF  SHIFTS 
!  GAUSSJ  RETURNS  THE  SHIFTS  AS  VECTOR  IN  DXIDV 
CALL  GAUSSJ(D2XIDADB,N  VAR,  1 00, DXIDV,  1 , 1 00) 

!  STORE  SHIFT  IN  DELA;  RETURNED  IN  DXIDV 
DO  I=1,NVAR 
DELA(I)=DXIDV(I) 

ENDDO 

!  UPDATE  THE  VARIABLES;  *******************APPLY  SHIFTS****************** 

!  WRITE(6,*)'VARIABLES  SHIFTS:' 

!  WRITE(  1 0,*)'VARI  ABLES  SHIFTS 

DO  I=1,NVAR 

!  WRITE(6, 303)1, V(I),DELA(I) 

!  WRITE(  10, 303)1,  V(I),DELA(I) 

303  FORMAT(lX,'#',I2,'  VAR=',F9.5,'  SHIFT=',F9.5) 

V(I)=V(I)+DELA(I) 

!  DAMP  KINETIC  RATEK'S 
IF(LLE.16)GOTO  304 
V(I)=V(I)-DELA(I)+0.60D0*(DELA(I)) 

304  CALL  VARFIG(I,NCOM,TAU,RATEK,V) 

ENDDO 

!  ESTIMATE  RESTRAINTS  BASED  ON  ALOW,AHI 
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DO  K=1,IMPTIVEC 
BLO(K)=l  .OD+25 
BHI(K)=1  .OD-25 
ENDDO 
DO  1=1, IR 
DO  K=1,IMPTIVEC 
SP=0.0D0 

DO  L=1,IMPTIVEC 
SP=SP+TAU(K,L)*ROW(I,L) 

ENDDO 

IF(SP.LT.BLO(K))BLO(K)=SP 

IF(SP.GT.BHI(K))BHI(K)=SP 

ENDDO 

ENDDO 

DO  I=1,IMPTIVEC 
DELOB  S= AHI-ALOW 
DELCALC=BHI(I)-BLO(I) 

SCALE(I)=DELCALC/DELOBS 
WRITE(10, 387)1, BLO(I),BHI(I), SCALE® 

387  F0RMAT(1X, 'SPECTRUM ',13,' LOW ',F8.3,' TO ',F8.3,';  SCALE=’,F8.3) 
ENDDO 

!  APPLY  RESTRAINTS  AFTER  THE  SHIFTS  *IF  OUT  OF  RANGE* 

DO  L=1,IMPTIVEC 
IF(SCALE(L).LE.1.3DO)GOTO  391 
SCALE(L)=(1.0D0-NCY/ICY)+((NCY/ICY)*SCALE(L)) 

DO  K=1,IMPTIVEC 
TAU(K,L)=TAU(K,L)/SCALE(L) 

ENDDO 
391  ENDDO 

!  GO  BACK  AND  ITERATE 
GOTO  285 

!  END  LOOP  OVER  ALL  VARIABLES 
!  END  OF  BIG  VARIABLE  ITERATION  LOOP 
!  NV  VARIABLES,  NK  RATE  CONSTANTS 
415  NT=NVAR-NK 
NT1-NT+1 
WRITE(  10,422) 

WRITE(  1 0,424)(V(I),I=1  ,NT) 

WRITE(  10,423) 

WRITE(  1 0,424)(V(I),I=NT  1  ,N  VAR) 

WRITE(6,422) 

WRITE(6,424)(V(I),I=1  ,NT) 

WRITE(6,423) 

WRITE(6,424)(V(I),I=NT1  ,NVAR) 

422  F0RMAT(1X, 'VARIABLES:  TRANSFORM’) 

423  F0RMAT(1X, 'VARIABLES:  RATE  CONSTANTS') 

424  FORMAT  ( 1  X,4(E  1 5 . 5 , 1 X)) 

!  SAVE  THE  CURRENT  VARIABLE  VALUES 
IF(IREAD.EQ.O)GOTO  994 
!  BACKSPACE  THE  POINTER 
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NR=0 

NLINES=(NVAR/3) 

FLINES=(NVAR/3.0D0) 

RESID=FLINES-NLINES 
IF(RESID.NE.0.0D0)NR=1 
NLIN  ES=NLINES+ 1 +NR 
DO  I=1,NLINES 
BACKSPACE(8) 

ENDDO 

994  WRITE(8,992)FUDGE 
WRITE(8,993)(V(I),I=1  ,NVAR) 

992  FORMAT(E15.6) 

993  FORMAT(3E15.6) 

WRITE(  1 0,*)'DATA  MATRIX  FIT  -  CALCULATED’ 

DO  1=1, IR 

WRITE(10,1 1)(ACALC(I,J),J=1,IC) 

ENDDO 
SSE=0.0D0 
DO  1=1, IR 
DO  J=1,IC 

!  Q(I,J)=Q(I,J)*(SIG(I,J)**2) 

SSE=SSE+Q(I,J)**2 

ENDDO 

ENDDO 

WRITE(  1 0,4 1 9)SSE 

419  FORM AT( IX, 'DIFFERENCE  MATRIX  SSE=’,E12.6) 
DO  1=1  ,IR 

WRITE(  1 0,428)(Q(1,J),J=1  ,IC) 

ENDDO 

428  FORMAT(1X,10(F6.4,1X)) 

!  CALCULATE  THE  ROTATED  SPECTRA 
999  CONTINUE 

WRITE(6,*)(TITLE(I),I=1 ,80) 

WRITE(  1 0,*)(TITLE(I),I=1 ,80) 

DO  I=l,NCOM 

WRITE(  10,861  )(CONC(I,J),J=  1  ,IC) 

861  FORMAT(1X,10(F6.4,1X)) 

ENDDO 

!  CALCULATE  THE  SPECTRA  (ROTATED  ROW  VECTORS) 
DO  1=1, IR 
DO  K=1,IMPTIVEC 
SP=0.0D0 

DO  L=1,IMPTIVEC 
SP=SP+TAU(K,L)*ROW(I,L) 

ENDDO 
G(I,K)=SP 
ENDDO 
ENDDO 
DO  1=1, IR 

WRITE(  1 0,865)WAV(I),(G(I,K),K=1  ,IMPTIVEC) 
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ENDDO 

865  FORMAT(1X,F4.0,1X,9(F7.4,1X)) 

1001  WRITE(6,*)"INSPECT  OUTPUT  ON  FILE  TESTVEC.OUT" 

CLOSE(8) 

CLOSE(IO) 

STOP 

END 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!  SUBROUTINES 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!  COMPUTE  THE  TRANSPOSE  OF  A  RECTANGULAR  MATRIX 
SUBROUTINE  TRANSPOSE, JJ, A, B) 

!  II  ROWS  AND  JJ  COLUMNS  OF  ORIGINAL  MATRIX:  JJ  ROWS,  II  COLUMNS  TRANSPOSED 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00, 1 00), B(  1 00, 1 00) 

!  WRITE(6,*)"TRANSPOSE",II,JJ 
DO  1=1,11 
DO  J=1,JJ 
B(J,I)=A(I,J) 

!  WRITE(10,*)B(J,I) 

ENDDO 

ENDDO 

RETURN 

END 


!  COMPUTE  THE  MATRIX  PRODUCT  OF  AN  IxK  AND  KxJ  ARRAYS  TO  PRODUCE  A 
!  THE  GENERAL  PRODUCT  MATRIX  OF  DIMENSION  IxJ 
SUBROUTINE  MATMULT(II,KK,JJ, A, B,AB) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00, 1 00), B(  1 00, 1 00),  AB( 1 00, 1 00) 

!  BY  THE  ROWxCOLUMN  METHOD;  AN  II  BY  KK  MATRIX  IS  MULTIPLIED  BY  A 
!  KK  BY  JJ  MATRIX  PRODUCING  AN  II  BY  JJ  MATRIX 
DO  1=1,11 
DO  J=1,JJ 
SUM=0.0D0 
DO  K=1,KK 

SUM=SUM+(A(I,K)*B(K,J)) 

!  WRITE(  1 0,40)1,  J,K,B(I,K),A(K,  J), SUM 

!  40  F0RMAT(1X,3I5,1X,3(F12.4,1X)) 

ENDDO 

AB(I,J)=SUM 

ENDDO 

ENDDO 

RETURN 

END 

] 

!  MULTIPLY  A  VECTOR  (B)  BY  A  MATRIX  (A)  [DOT  PRODUCT] 
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!  COLUMNS  OF  VECTOR  MUST  EQUAL  ROWS  OF  MATRIX  (IR) 
!  OUTPUT  VECTOR  (C)  HAS  IR  ROWS 
SUBROUTINE  MATVEC(IR,IC,A,B,C) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00, 1 00), B(  1 00), C(  1 00) 

DO  1=1, IC 
SUM=0.0D0 
DO  J=1,IR 

SUM=SUM+(A(I,J)*B(J)) 

!  WRITE(  1 0,30)I,J,A(I,J),B(J),SUM 
!  30  F0RMAT(1X,I5,I5,3(F12.5,1X)) 

ENDDO 

C(I)=SUM 

ENDDO 

RETURN 

END 

i 

!  MULTIPLY  TWO  VECOTRS  (A)  AND  (B);  (RxC)*(CxR)  =  (RxR) 

!  TO  PRODUCE  A  MATRIX 

SUBROUTINE  VECMULT(II,JJ,A,B,C) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(  1 00), B(  1 00),C(1 00, 1 00) 

DO  1=1,11 
DO  J=1,II 
C(I,J)=A(I)*B(J) 

!  WRITE(1 0,30)1, J,A(I),B(J),C(I,J) 

!  30  F0RMAT(1X,2I5,1X,3F12.5) 

ENDDO 

ENDDO 

RETURN 

END 


*  *  *  *  *  *  *  sfc  *  *  *  * 

MINV 

SUBROUTINE  MINV(AB,N,ND, SCRATCH, DET, EPS, M, MODE) 
IMPLICIT  REAL*8(A-H,0-Z) 


!  A  subroutine  that  calculates  the  determinant  and  inverse  of 
!  a  matrix,  as  well  as  solving  systems  of  linear  equations. 

!  Martin  J.  McBride.  11/25/85. 

!  General  Electric  CRD,  Information  System  Operation. 

INTEGER  N,ND,M, MODE, OUTER, ROW, COL, I, SCOL,SROW,PIVCNT 
DIMENSION  AB(  1 00, 1 00),SCRATCH(500) 

!  Initialize  scratch  space,  with  1  to  N  holding  the  diagonal  of  the  identity 
!  matrix  used  to  compute  the  inverse  and  N+l  to  2N  holding  the  positions  of 
!  the  first  N  columns  of  the  matrix  (for  use  when  pivot  occurs). 

DO  5  I  =  1  ,N 
5  SCRATCH(I)  =  1.0 
COLNUM=  1.0 
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DO  6  I  =  N+1,2*N 
SCRATCH(I)  =  COLNUM 
COLNUM  =  COLNUM  +  1 .0 
6  CONTINUE 

!  Make  left,  square  matrix  an  upper  triangular  matrix. 

DET  =  0.0 

PIVCNT  =  0 

DO  10  OUTER  =  1,N-1 

IF  (DABS(AB(OUTER, OUTER))  .LE.  EPS)  THEN 
CALL  PIVOT(AB,N,ND, OUTER, SCRATCH, EPS) 

IF  (AB(OUTER, OUTER)  .EQ.  0.0)  THEN 
WRITE(6,*)'*************************************t 
WRITE(6,*)'  MINV  called  with  singular  matrix.' 
WRITE(6,*)'*************************************i 
STOP 
ENDIF 

PIVCNT  =  PIVCNT  +  1 
ENDIF 

DO  20  ROW  =  OUTER+l,N 
MULT  =  AB(ROW,OUTER)/AB(OUTER, OUTER) 

DO  30  COL  =  OUTER, N+M 

30  AB(ROW,COL)  =  AB(ROW,COL)  -  AB(OUTER,COL)*MULT 
DO  25  SCOL  =  1 , OUTER- 1 

25  AB(ROW,SCOL)  =  AB(ROW,SCOL)  -  AB(OUTER,SCOL)*MULT 
AB(ROW, OUTER)  =  AB(ROW, OUTER)  -  SCRATCH(OUTER)*MULT 
20  CONTINUE 
10  CONTINUE 
!  Compute  determinant. 

DET  =  AB(1,1) 

DO  40  I  =  2,N 
40  DET  -  DET*AB(I,I) 

DET  =  (- 1 .0)*  *PI  VCNT  *  DET 

!  Return  if  inverse  is  not  to  be  found  and  there  are  no  systems  of  equations 
!  to  solve. 

IF  (MODE  .EQ.  0  .AND.  M  .EQ.  0)  RETURN 
!  Place  ones  in  diagonal  of  square  matrix  A. 

DO  80  ROW  =  1  ,N 
DIV  =  AB(ROW,ROW) 

DO  90  COL  =  1,N+M 
AB(ROW,COL)  =  AB(ROW,COL)/DIV 
90  CONTINUE 

SCRATCH(ROW)  =  SCRATCH(ROW)/DIV 
80  CONTINUE 

•'  Reduce  upper  triangle  to  zeros  to  give  matrix  A  =  I 
DO  50  OUTER  =  2,N 
DO  60  ROW  =  OUTER- 1, 1,-1 
MULT  =  AB(ROW,OUTER)/AB(OUTER,  OUTER) 

DO  70  COL  =  OUTER,N+M 

70  AB(ROW,COL)  =  AB(ROW,COL)  -  AB(OUTER,COL)*MULT 
DO  65  SCOL  =  1 , ROW-1 
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65  AB(ROW,SCOL)  =  AB(ROW,SCOL)  -  AB(OUTER,SCOL)*MULT 
SCRATCH(ROW)  =  SCRATCH(ROW)  -  AB(OUTER,ROW)*MULT 
DO  63  SCOL  =  ROW+1, OUTER- 1 

63  AB(ROW,SCOL)  =  AB(ROW,SCOL)  -  AB(OUTER,SCOL)*MULT 

AB(ROW, OUTER)  =  AB(ROW, OUTER)  -  SCRATCH(OUTER)*MULT 
60  CONTINUE 
50  CONTINUE 

!  Move  diagonals  of  inverse  to  matrix  AB. 

DO  85  I  =  1,N 

85  AB(I,I)  =  SCRATCH(I) 

!  If  pivot  was  made,  switch  rows  corresponding  to  the  columns  that  were 
!  pivoted. 

IF  (PIVCNT  .EQ.  0)  RETURN 
ROW  =  1 
DO  95  I  =  1,N-1 

SROW  =  INT(SCRATCH(ROW+N)) 

IF  (SROW  .NE.  ROW)  THEN 
DO  92  COL  =  1  ,N+M 
TEMP  =  AB(ROW,COL) 

AB(ROW,COL)  =  AB(SROW  ,COL) 

AB(SROW,COL)  =  TEMP 
92  CONTINUE 

TEMP  =  SCRATCH(ROW+N) 

SCRATCH(ROW+N)  =  SCRATCH(SROW+N) 

SCRATCH(SROW+N)  =  TEMP 
ELSE 

ROW  =  ROW  +  1 
ENDIF 

95  CONTINUE 
RETURN 
END 


SUBROUTINE  PIVOT(AB,N,ND, OUTER, SCRATCH, EPS) 

IMPLICIT  REAL*8(A-H,0-Z) 

This  subroutine  switches  two  columns  of  a  matrix  to  get 
a  nonzero  entry  in  the  diagonal. 

Martin  J.  McBride.  12/04/85. 

General  Electric  CRD,  Information  System  Operation. 

INTEGER  N  ,ND, COL, OUTER, I 
DIMENSION  AB(100,100),SCRATCH(100) 

!  Get  first  column  with  non-zero  element  in  row  OUTER. 

COL  =  OUTER  +  1 
10  IF  (COL  .GT.  N)  GO  TO  90 
IF  (ABS(AB(OUTER,COL))  .GT.  EPS)  GO  TO  20 
COL  =  COL  +  1 
GOTO  10 

i  Switch  column  OUTER  with  column  COL,  which  has  non-zero  element  in 
1  row  OUTER. 
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20  DO  301  =  1,N 

TEMP  =  AB(I, OUTER) 

AB(I,  OUTER)  =  AB(I,COL) 

AB(I,COL)  =  TEMP 
30  CONTINUE 

TEMP  =  SCRATCH(N+OUTER) 

SCRATCH(N+OUTER)  =  SCRATCH(N+COL) 

SCRATCH(N+COL)  =  TEMP 
90  CONTINUE 
RETURN 
END 

! 

SUBROUTINE  ALSQ(N,X,Y,  SLOPE,  ERSLOP,YINT,ERYINT,CORREL,SSE) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  X(100),Y(100) 

NT=N-2 

IF(NT.EQ.O)GO  TO  80 

SUMX=0.0D0 

SUMY=0.0D0 

SUMXX=0.0D0 

SUMXY=0.0D0 

SUMYY=0.0D0 

SSE=0.0D0 

DO  1=1, N 

SUMX=SUMX+X(I) 

SUM  Y=SUMY +Y  (I) 

SUMXX=SUMXX+(X(I)  *X(I)) 

SUM  Y  Y=SUM  Y  Y +( Y  (I)  *  Y(I)) 

SUMXY=SUMXY+(X(I)*Y  (I)) 

ENDDO 

SLOPE=((N*  SUMXY)-(SUMX*  SUMY))/ ((N*SUMXX)-SUMX*  *2) 

!  WRITE(6,*)SLOPE 

YINT=(SUMY/N)-(SLOPE*(SUMX/N)) 

DO  1=1  ,N 

SSE=SSE+((SLOPE*X(I))+YINT-Y(I))**2 

ENDDO 

DEL=(N*  SUMXX)-(SUMX*  SUMX) 

SIG2=(1.0D0/(N-2.0D0))*SSE 

ERSLOP=DSQRT((N*SIG2)/DEL) 

ERYINT=DSQRT(SIG2*SUMXX/DEL) 

CNUM=(N*  SUMXY)-(SUMX*  SUM  Y) 
TERM1=DSQRT((N*SUMXX)-SUMX**2) 

TERM2=DSQRT  ((N*  SUMYY)-SUMY  *  *2) 

CORREL=CNUM/(TERMl  *TERM2) 

!  WRITE(6,*)CORREL 
GO  TO  90 
80  WRITE(6,82) 

82  FORMAT(  IX, 'INSUFFICIENT  DATA  FOR  REGRESSION') 

90  RETURN 
END 
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SUBROUTINE  VARFIG(L,NCOM,TAU,RATEK,V) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  TAU(  1 0, 1 0),RATEK(  1 0), V(  1 00) 

!  NUMBER  OF  VARIABLES  IS  NCOMA2  +  (NCOM-1) 

!  AFTER  DETERMINING  WHAT  VARIABLE  "L"  TRANSLATES  INTO,  IT  CORRECTS  ONLY 
THIS 

NA=NCOM 

NB=NCOM-l 

IDIV=(NCOM**2)+l 

IF(L.GE.IDIV)GOTO  396 

IX=(L+NB)/NA 

IY=L+1-((NA*IX)-NB) 

GOTO  397 

396  IZ=L-(IDIV-1) 

RATEK(IZ)=V  (L) 

GOTO  398 

397  TAU(IX,IY)=V(L) 

398  RETURN 
END 

! 

SUBROUTINE  CONCEN(NC,J,RATEK,TIME,CN,CONC) 

!  THIS  ROUTINE  FOR  A=B=C=D  (FOR  1ST  1ST  2ND  ORDER,  RESP)  7/1/00  EJV 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  RATEK(  1 0),TIME(  1 00),CN(  1 0),CONC(  1 0, 1 00) 

!  DIMENSION  CNUM(10),CDEN(  10) 

!  FIRST  SPECIES 

CN(  1  )=DEXP(-RATEK(  1  )*TIME(J)) 

!  SECOND  SPECIES 

B=(DEXP(-RATEK(  1  )*TIME(J)))-(DEXP(-RATEK(2)*TIME(J))) 

CN(2)=B*(RATEK(  1  )/(RATEK(2)-RATEK(  1 ))) 

!  THIRD  SPECIES 

CALL  RUNGE(J,CC,RATEK,TIME) 

CN(3)=CC 

!  AND  ADD  THE  NCth  SPECIES 
CAL=0.0D0 
DO  1=1,3 
C  AL=C  AL+CN  (I) 

ENDDO 

CN(4)=1 .0D0-CAL 

! 

DO  1=1,  NC 
CONC(I,J)=CN(I) 

ENDDO 

RETURN 

END 

! 

SUBROUTINE  RUNGE(J,CC,RATEK,TIME) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  RATEK(  1 0),TIME(  1 00) 
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!  DIMENSION  DELC(100),DCDT(100) 

A 1  =(RATEK(  1)*RATEK(2))/(RATEK(2)-RATEK(  1 )) 

!  SPECIFY  TIME  INTERVAL  (TINT) 

TINT=10.0D0 
IT=TIME(J)/TINT 
TIM=0.0D0 
CONC=O.ODO 
DO  1=1,  IT 

!  GENERATE  THE  FOUR  RUNGE-KUTTA  TERMS 
CALL  KUTTA(TIM,CONC,  A 1  ,TINT,RATEK,  A) 

RK1=A 

TH=TIM+TINT/2.0D0 
CH=CON  C+A/2 . 0D0 

CALL  KUTTA(TH,CH,A1,TINT,RATEK,B) 

RK2=B 

CI=CONC+B/2.0D0 

CALL  KUTTA(TH,CI,A1  ,TINT,RATEK,C) 

RK3=C 

TI=TIM+TINT 

CJ=CONC+C 

CALL  KUTTA(TI,CJ,A1,TINT,RATEK,D) 

RK4=D 

DELC=RK1  +(RK2*2 ,0D0)+(RK3  *2.0D0)+RK4 

CONC=CONC+DELC 

TIM=TIM+TINT 

ENDDO 

CC=CONC 

RETURN 

END 

t 

SUBROUTINE  KUTTA(T,C,A1  ,TINT,RATEK,A) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  RATEK(IO) 

A=(A1*DEXP(4RATEK(1)*T)-A1*DEXP(-RATEK(2)*T)-RATEK(3)*(C**2))*TINT 

RETURN 

END 

i 

SUBROUTINE  MODEL(I,NCOM,CN,TAU,ZETA, ROW, TERM) 

!  THIS  ROUTINE  MADE  GENERAL  FOR  UP  TO  9  NCOM  COMPONENTS  6/24/00  EJV 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  CN(  1 0),TAU(  10,1 0),ZETA(  10,1 0),ROW(  1 00, 1 00) 

DO  L=l,NCOM 
DO  K=l,NCOM 
ZETA(K,L)=TAU(L,K)*CN(L) 

ENDDO 

ENDDO 

!  ACCUMULATE  THE  TERMS  IN  ABSORBANCE 
TERM=0.0D0 
DO  L=l,NCOM 
DO  K=l,NCOM 
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TERM=TERM+(ROW(I,K)*ZETA(K,L)) 

ENDDO 

ENDDO 

RETURN 

END 

! 

SUBROUTINE  GAUSSJ(A,N,NP,B,M,MP) 

!  THIS  ROUTINE  SOLVES  A  x  =  B  BY  (AT.A).x=(AT.B)  (AT=Atranspose) 

!  IT  TAKES  MATRIX  A(N,N)  STORED  IN  AN  ARRAY  OF  DIMENSIONS  NPxNP 
!  B  IS  AN  INPUT  MATRIX  NxM  CONTAINING  m  RIGHT-SIDE  VECTORS,  STORED  IN  AN 
!  ARRAY  OF  DIMENSIONS  NPxMP.  ON  OUTPUT,  A  IS  REPLACED  BY  ITS  MATRIX 
!  INVERSE,  AND  B  IS  REPLACED  BY  THE  CORRESPONDING  SET  OF  SOLUTION  VECTORS 
IMPLICIT  REAL*8(A-H,0-Z) 

!  PARAMETER  (NMAX=50) 

DIMENSION  A(  1 00, 1 00), B(  1 00, 1 00),IPIV(50),INDXR(50),INDXC(50) 

DO  11  J=1,N 
IPIV(J)=0.0D0 

11  CONTINUE 
DO  22  1=1, N 
BIG=0.0D0 
DO  13  J=1,N 
IF(IPIV(J).NE.1)THEN 
DO  12  K=1,N 
IF(IPIV(K).EQ.0)THEN 
IF(DABS(A(J,K)).GE.BIG)THEN 
BIG=DABS(A(J,K)) 

IROW=J 

ICOL=K 

ENDIF 

ELSE  IF  (IPIV(K).GT.l)  THEN 
PAUSE  'SINGULAR  MATRIX' 

ENDIF 

12  CONTINUE 
ENDIF 

13  CONTINUE 

IPI  V(ICOL)=IPI  V(ICOL)+ 1 
IF(IROW.NE.ICOL)  THEN 
DO  14  L=1,N 
DUM=A(IROW,L) 

A(IROW,L)=A(ICOL,L) 

A(ICOL,L)=DUM 

14  CONTINUE 
DO  15  L=1,M 
DUM=B(IROW,L) 

B(IROW,L)=B(ICOL,L) 

B(ICOL,L)=DUM 

15  CONTINUE 
ENDIF 

INDXR(I)=IROW 

INDXC(I)=ICOL 
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IF ( A(ICOL,ICOL)  .EQ . 0. ODO)  PAUSE  'SINGULAR  MATRIX' 
PIVINV=1  ,ODO/A(ICOL,ICOL) 

A(ICOL,ICOL)=l  .ODO 
DO  16,L=1,N 

A(ICOL,L)=A(ICOL,L)*PIVINV 

16  CONTINUE 
DO  17  L=1,M 

B(ICOL,L)=B(ICOL,L)*PIVINV 

17  CONTINUE 
DO  21  LL=1,N 
IF(LL.NE.ICOL)THEN 
DUM=A(LL,ICOL) 

A(LL,ICOL)=O.ODO 
DO  18  L=1,N 

A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 

18  CONTINUE 
DO  19  L=1,M 

B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 

19  CONTINUE 
ENDIF 

21  CONTINUE 

22  CONTINUE 
DO  24  L=N,1,-1 

IF(INDXR(L).NE.INDXC(L))THEN 
DO  23  K=1,N 
DUM=A(K,INDXR(L)) 

A(K,IN  DXR(L))=A(K,IN  DXC(L)) 

A(K,INDXC(L))=DUM 

23  CONTINUE 
ENDIF 

24  CONTINUE 
RETURN 
END 
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First  Order  Kinetics  Experiments 

Global  analysis . Version  1.1  EJV  08/09/2000 

Principal  component  analysis  (PCA)  method  -  HP8453  Array  Detector  Data 


Obsrve  the  following  data  analysis  scheme.  The  FORTRAN  programs  are  on  chelablO  under 
C:/LF90/CR/.  The  programs  are  preparation  for  principal  component  analysis  (PREPPCA),  principal 
component  analysis  (PCA),  and  vector  testing  and  least-squares  (TESTVEC). 

Scheme:  activities,  PROGRAMS,  files 
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Detailed  instructions. 

Activity. 

HP8453  Diode  Array  Spectometer.  Kinetic  Run. 

1 .  Set  up  a  kinetics  run.  Use  the  jacketted  cell,  a  thermostated  water  bath  with  recirculating  pump.  Give 
the  kinetics  run  datafile  an  appropriate  name. 

2.  For  a  25  °C  run,  use  a  25-s  ‘start’  time,  and  collect  spectra  every  60  s  thereafter  for  a  period  of  3,600  s. 
Adjust  this  for  higher  or  lower  temperatures. 

3.  Combine  thermal-equilibrated  reagents  at  0  time  and  initiate  the  kinetics  run;  mix!;  transfer  to  the  cell 
within  25  s.  Sit  back  and  collect  data.  Whole  spectra  are  stored  in  the  kinetics  run  file. 

Data  processing  and  transfer. 

4.  Selecting  the  individual  spectra  with  the  mouse  (right  click  on  spectra),  save  the  spectra  to  be  analyzed 
as  .dif  files  to  diskette.  Note  the  times  of  the  spectra.  Each  file  should  be  saved  as  a  filename  that  makes 
sense  with  regard  to  the  time,  like  A85.dif  for  the  file  containing  the  data  from  the  85-s  spectrum. 

5.  After  selecting  and  saving  as  .dif  files  all  the  spectra  of  interest,  place  them  in  C:/LF90/CR7  on 
chelablO  (MSDOS  prompt). 

6.  Run  PREPPCA  (invoke  PREPPCA)  to  prepare  the  data  for  principal  component  analysis.  In  this 
program,  you  will  be  asked  to: 

a)  Enter  a  range  of  wavelengths  (in  nm)  to  be  analyzed  (like  400,  700  [return]) 

b)  Enter  the  uniform  interval  at  which  the  spectra  are  to  be  tabulated  (like  10  [return]) 

c)  Enter  the  name  of  the  .dif  files  in  order  of  the  times  collected  (like  A85.dif  [return])  followed  by  the 
time  (in  seconds)  (like  85.  [return]);  and  repeat  this  until  all  of  the  data  have  been  entered.  Then  enter 
STOP  [return].  If  you  make  an  error,  execute  <control>C,  then  rerun  the  program. 

d)  Enter  a  name  for  the  data  file  (like  Cr+6  to  Cr+3,  temperature,  date,  name,  class  [return]). 

7.  PREPPCA  runs  and  outputs  three  files:  hptrial.dat,  preppca.out,  and  info.out. 

Inspect  the  hptrial.dat  file  to  see  that  it  contains: 

a)  The  data  matrix  preceded  by  the  number  of  rows  (wavelengths)  and  columns  (times). 

b)  The  wavelengths  and  then  the  times.  Inspect  these  carefully  to  be  sure  they  are  correct. 

c)  The  title  of  the  kinetics  run. 

If  all  is  well,  change  the  name  of  this  file  to  hpdata.dat.  This  will  be  an  input  file  to  PCA  (principal 
component  analysis).  The  file  preppca.out  contains  similar  information,  and  can  be  printed  out  if  desired. 

8.  The  file  info.out  contains  representations  of  the  first  and  last  spectra  which  will  be  useful  for  vector 
testing  later  on.  In  fact,  for  a  two-component  system,  the  first  spectra  should  be  the  Cr+6  spectrum,  the 
second  should  be  the  Cr+3  spectrum.  The  last  is  approximated  by  the  t  =  large  spectrum  collected;  the 
Cr+6  spectrum  is  estimated  from  the  first  spectrum  back  extrapolated  to  t  =  0. 

Principal  component  analysis. 

9.  There  is  no  program  control  for  principal  component  analysis.  Be  sure  that  its  input  data  are  present  in 
file  hpdata.dat.  Execute  the  program  by  invoking  PCA. 

10.  PCA  outputs  two  files:  pca.out  and  testvec.out.  The  file  pca.out  can  be  inspected  or  printed  out 
(partially,  since  it  is  possibly  a  large  file)  to  answer  the  following  questions: 

a)  How  many  factors  are  needed  to  explain  the  variance  in  the  kinetics  data? 

b)  How  large  are  the  eigenvalues,  and  how  many  are  there  before  their  magnitude  vanishes? 

c)  Do  the  row  vectors  resemble  the  spectra  of  the  components? 

d)  Do  the  column  vectors  resemble  the  concentration  dependence  of  the  components  over  time? 

If  the  data  can  be  explained  with  two  components,  then  proceed  to  test  spectral  vectors  and  resolve  the 
kinetics  information. 
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Vector  testing  and  least  squares  analysis. 

1 1 .  Examine  the  info. out  file.  At  the  end  of  the  file,  control  information  for  TESTVEC6  is  required.  It  is 
contained  in  columns  1  -  5  on  the  line  immediately  after  the  second  test  vector. 

These  are  ILSQ,  NCOM,ICONT,IREAD,IFINISH  (in  511  format;  see  below). 

Recommended  values  at  the  start  of  a  run  are:  00000 
Rename  the  file  to  info.dat. 

ILSQ  =  0  do  not  perform  least  squares;  =  1  perform  least  squares. 

NCOM  =  number  of  components;  default  (0)  is  the  number  of  important  vectors  (or  2). 

ICONT  =  run  2A(NCOM)  cycles  of  least  squares  (if  ILSQ  =  1);  default  NCOM  =  0  (i.e.,  1  cycle). 

IREAD  =  0  do  not  read  previous  info.dat  file  results;  =  1  read  previous  info.dat  file  results  and  begin  least 
squares  calculation  from  there,  and  write  over  old  info.dat  file  after  cycling. 

IFINISH  =  0  do  not  output  spectral  information  at  a  particular  wavelength;  =  1  prompt  at  end  of  the  run 
for  a  wavelength  (tanax)  at  which  to  output  time  vs.  absorbance  data. 

12.  Execute  TESTVEC6  (using  00000).  This  program  evaluates  the  test  vectors  contained  on  info.dat  file. 
It  output  is  in  one  file:  testvec.out. 

13.  Examine  the  file  testvec.out.  Answer  the  following  questions: 

a)  Are  the  two  test  vectors  likely  to  be  real  spectral  vectors?  (What  are  the  sum-of-square-errors  for  the 
test  and  least  squares  trial?) 

b)  Do  the  transforms  of  the  column  vectors  look  like  improved  concentration  vs.  time  data? 

14.  Now,  fit  the  transformation  matrix  elements  and  a  single  rate  constant  by  the  nonlinear  least  squares 
method.  This  is  done  by  reinvoking  TESTVEC6  with  the  info.dat  command  of  lOnlO  (where  n  =  2An 
cycles;  10410  will  run  16  cycles  and  save  the  results).  Rerun  until  converged. 

15.  Examine  the  file  testvec.out.  Answer  the  following  additional  questions: 

a)  What  rate  constant  accounts  for  the  data? 

b)  How  well  does  the  kinetics  model  reproduce  the  spectral  data?  Are  the  residuals  small? 

16.  Now,  end  the  process  by  outputing  the  spectal  information  at  the  kmax  (say,  440nm).  This  is  done  by 
reinvoking  TESTVEC6  with  the  info.dat  command  of  101 1 1.  A  prompt  at  the  end  of  the  program  will 
inquire  for  a  wavelength,  and  the  output  will  be  on  file  testvec.out.  This  can  be  clipped  out  and  plotted 
with  the  other  spectral  and  time  data. 
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